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ATTITUDE  CQMTROL  SYSTEM  OF  RECOVERABLE 
EARTH-ORIENTING  TECHNOLOGY  TESTING  SATELLITES 
AND  SOME  FLIGHT  TESTING  RESULTS 


Yang  Jiachi,  Zhang  Guofu,  Sun  Chengqi,  Feng  Xueyi 
and  Niu  Yinsheng* 

(Beijing  Institute  of  Control  Engineering) 

ABSTRACT 

In  1970,  three  recoverable  technology  testing  satellites, 
permitting  orientation  of  observatory  with  respect  to  the 
earth,  were  successfully  launched,  ooeratea  in  orbit  and 
recovered.  This  paper  reviews  the  three-axis  stabilized 
attitude  control  system  used  in  these  satellites,  including 
its  composition,  function  and  considerations  of  selecting 
main  parameters  in  the  system,  etc.  Finally  the  opening 
performances  of  the  flight  testing  results  are  briefly  described. 

INTRODUCTION 


A  first  generation  three-axis  stabilized  attitude  control 
system  was  first  used  on  recoverable  technology  testing  satellites 
permitting  orientation  of  observation  with  respect  to  the  Earth 
and  successfully  accomplished  the  following  tasks: 

Eliminating  the  initial  satellite  attitude  aberration  caused  by 
the  separation  of  the  launch  vehicle  from  the  satellite: 

assuring  the  attitude  precision  and  the  attitude  angular  velocity 
reauirement  needed  for  satellite  observation  of  the  Earth  during  the 
crbitting  period; 

enabling  the  satellite  to  achieve  the  attitude  required  for 
recovery  by  rapidly  rotating  about  the  yaw  axis  before  returning  to 
Earth . 

In  this  paper,  we  emphasize  the  introduction  of  the  structure, 
the  mechanism  and  the  design  for  the  principal  parameters  of  the 
attitude  control  system.  Finally,  the  analysis  of  the  result  of  the 
flight  tests  is  briefly  described. 
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II.  DESCRIPTION  OF  THE  ATTITUDE  CONTROL  SYSTEM  AMD  THE  CONSIDERATION 
FOR  THE  CHOICE  OF  THE  PRINCIPAL  PARAMETERS 

The  attitude  control  system  is  formed  of  the  attitude  measure¬ 
ment  system,  the  central  control  circuits  and  the  propulsion  execution 
mechanism  through  the  satellite  closed  circuit.  The  description  of 
the  attitude  movement  is  based  upon  the  Earth-centered  orbital  coordin¬ 
ate  system  and  the  satellite  coordinate  system  as  shown  in  Figure  1. 

The  angles  between  the  two  systems  are  denoted  as  t,  fl  and  respec¬ 
tively  called  the  roll  attitude  angle,  the  pitch  attitude  angle  and 
the  yaw  attitude  angle. 

1.  The  attitude  measurement  system 

In  order  to  establish  the  Earth-centered  orbital  coordinate 
system  in  the  orbit  as  a  standard  for  measuring  the  satellite  attitude 
and  to  give  the  signal  for  attitude  deviation  while  being  able  to 
satisfy  the  requirement  of  eliminating  the  initial  attitude  devia¬ 
tion  and  establishing  the  return  attitude,  we  have  selected  an  atti¬ 
tude  measurement  system  composed  of  two  conical  scanning  infrared 
level  meters  and  two  gyros  with  two-three  decrees  of  freedom.  As 

shown  in  Figure  2,  seme  improvement  to  the  characteristics  of  the 
level  meter  has  been  adopted  in  the  design  of  the  optical  system  and 
the  electronic  circuitry  of  the  level  meter.  The  complementary 
characteristics  of  the  infrared  level  meters  and  the  gyro  have  been 
fully  utilised  in  this  measurement  system.  The  sr.yrc  output  has  a  small 
high  frequency  noise  but  also  a  significant  constant  drift.  Hence, 
it  cannot  be  used  as  a  long  term  standard  of  measurement  by  itself 
while  the  infrared  level  meter  has  a  small  constant  drift.  Thus  it  is 
usable  as  a  long  term  standard  for  measurement,  but  it  has  significant 
alternating  noise.  It  cannot  give  the  measurement  standard  for  the 
yaw  altitude  and  cannot  be  used  by  itself.  By  combining  these  two 
elements,  they  will  complement  each  other  and  form  a  better  attitude 
measurement  system. 

From  Figure  2,  it  can  be  seen  that  the  long  term  measurement 
standard  of  this  measurement  system  is  provided  by  the  infrared 
level,  and  the  transient  measurement  standard  is  provided  by  the  gyro 
with  three  channels  of  relatively  smooth  attitude  deviation  signals 
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Figure  2.  Schematic  of  the  attitude  measurement  system  composition. 

• —  transducer  • —  force  coupling  device 

1 — pitch  angle  signal;  2 — pitch  infrared  levelmeter;  3 — amplification 
and  correct  ion ;  4 — phase  sensor;  5  —  compensation  signal;  6  —  large  angle 
transducer;  7 — pitch  axis  measuring  system;  8 — horizontal  gyro; 

9 — programmed  mechanism;  10 — roll  axis;  11 — pulse  source;  12 — pitch 
axis;  13 — amplification  and  correction;  14 — phase  sensor;  15 — yaw 
axis;  16 — yaw  angle  signal;  17 — phase  sensor;  18 — roll  and  yaw  axis 
measuring  system  (gyro  compass);  19 — vertical  gyro;  20 — amplification 
and  correction;  21 — roil  infrared  levelmeter;  22 — amplification  and 
correction;  23 — phase  sensor;  24 — roll  angle  signal 

3  a 


during  the  establishment  of  the  return  attitude-  The  large  angle 
transducer  may  be  used  to  monitor  the  working  condition  of  the  atti¬ 
tude  control  system  during  the  establishment  of  the  return  attituie- 
In  eliminating  the  initial  attitude  deviation  as  well  as  in  estab¬ 
lishing  the  return  attitude,  the  infrared  level  must  be  cut  off  from 
the  attitude  measurement  system  so  that  the  system  is  in  an  open 
circuit  measurement  mode,  i.e.,  the  two  gyros  are  both  situated  in  a 
free  gyro  state.  Figure  2  shows  the  working  status  of  the  measure¬ 
ment  system  in  orbit. 

When  the  orbit  is  nearly  circular.  Figure  2  may  be  simplified 
to  the  block  diagram  as  shown  in  Figure  3  within  the  linear  opera¬ 
tional  range  of  small  signals. 

From  Figure  3,  we  can  write  the  following  equations: 

i> )  =  K  *(4>h—  <f>c)  +  D, 


where  9,  \p  and  <J>  are  the  pitch,  yaw  and  roll  attitude  angles  of  the 
satellite.  0,,,  and  are  the  pitch,  yaw  and  roll  frame  angles 
of  the  gyro.  K0,  and  are  the  amplifications  on  the  gyro  advanc¬ 
ing  angular  velocities  provided  by  the  amplifier  and  the  gyro  coupling 
device.  D^,  and  are  the  equivalent  drift  velocities  of  the  gyro 
along  the  roll,  pitch  and  yaw  axes,  including  the  components  and 
quadratic  coupling  terms  along  the  corresponding  axes  of  the  gyro 
drift  velocity,  the  orbital  plane  advancing  angular  velocities,  etc., 
ft  and  ftQ  are  respectively  the  instantaneous  and  normalized  orbital 
angular  velocities,  the  latter  being  provided  by  the  amplifier.  Both 
Figure  3  and  equations  (1),  (2)  and  (3)  indicate  that  the  critical 
parameters  here  are  Kg,  and  K^.  The  principle  of  their  selection 
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Figure  3*  Simplified  block  diagram  of  the  attitude 

measuring  system. 

1 — pitch  axis  attitude  measuring  system;  2 — the  gyro  compass  system 

should  be  based  upon  the  noise  of  the  levelmeter  and  the  drift  rate 
of  the  gyro  on  the  corresponding  axes,  so  as  to  minimize  the  mean 
square  value  of  the attitucje  measurement  error.  At  the  same  time, 
consideration  must  be  given  to  the  effect  of  the  noise  of  the  level- 
meter  on  the  rate  of  gas  consumption  produced  through  the  central 
circuit.  This  is  important  to  satellites  of  medium  and  low  orbits, 
especially  when  infrared  levelmeter  wavelength  cannot  be  located  in 
the  Idealized  14-16  micron  band.  In  the  design  of  this  system,  we 
considered  the  maximum  drift  rate  of  the  gyro  is  4°  per  hour,  the 
mean  square  value  of  the  Infrared  levelmeter  error  is  0.25°,  1000 
times  the  bandwidth  for  the  steady  noise  for  the  orbital  angular 
velocity.  The  value  range  of  the  three  parameters  mentioned  above 
is  determined  under  these  conditions. 

The  following  problems  should  also  be  considered  in  the  system 
design: 
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(1)  Limiting  the  maximum  saturated  advancing  angular  velocity 
of  the  gyro:  To  avoid  the  increase  of  the  attitude  deviation  and 
excessive  system  gas  usage  rate  produced  by  the  high  output  signals 
when  the  levelmeter  scans  the  sun  or  the  moon  during  normal  flight 
of  the  satellite,  the  range  of  values  of  Kg,  and  should  be 
limited.  In  our  system  design,  we  require  that: 

the  saturated  pitch  advancing  angular  velocity  =  0.06% 
the  saturated  yaw  advancing  angular  velocity  ^  *  0.054% 
the  saturated  roll  advancing  angular  velocity  .<  =  0.018% 

(2)  The  problem  of  matching  the  transmission  coefficients  of 
the  infrared  levelmeter  and  the  gyro:  During  technical  implementa¬ 
tion,  the  transmission  coefficient  between  the  input  angle  and 
output  potential  of  the  infrared  levelmeter  may  deviate  from  the  trans¬ 
mission  coefficient  K  between  the  input  angle  of  the  gyro  and  the 

g 

transducer  output  potential.  Analysis  and  experiment  demonstrate  that 

when  Kh  is  different  from  K  ,  the  system  will  be  induced  to  operate 

in  a  limiting  loop  of  multiple  gas  propulsion,  thus  increasing  the  gas 

consumption  rate.  The  effect  on  the  roll  channel  is  more  pronounced 

than  the  other  two  channels.  For  example,  when  Kh  is  larger  than 

by  10%,  the  yaw  mini-propulsion  system  undergoes  a  steady  "double 

propulsion  phenomenon".  The  limiting  loop  angular  velocity  and  the 

gas  consumption  rate  both  double.  In  our  design,  we  require 

roll  channel  (K,  -  K  )/K  =  6% 

h  g  g 

yaw  and  pitch  channel  (K.  -  X  )/K  <  =  10% 

h  g  g 


(3)  Selection  of  the  gyro  phase  sensitivity  time  constant  T  . 

Apparently  the  increase  of  T  will  decrease  the  resistance  of  the 

S 

system.  When  T  reaches  a  certain  value,  the  joint  operation  of  the 
g 

large  and  small  propulsions  of  the  system  may  be  damaged.  The  dynamic 

characteristics  of  the  pseudovelocity  increment  feedback  system 

indicatethat  Tg  is  not  greater  than  its  recharge  time  constant  T^c . 

In  general,  we  should  let  T  *  =  0 . 2  . 

g 


(4) 


The  time  constant  T,  of  the  infrared  levelmeter: 

h 


The 


increase  of  T^  will  induce  the  steady  state  limiting  loop  operating 
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Figure  4.  Flow  chart  of  the  pitch  channel  control  system 
I  =  satellite  pitch  axis  rotational  inertia 
S  =  Laplacian 

1 — high  powered  execution  mechanism;  2 — low  powered  execution 
mechanism;  3 — attitude  measurement 

speed  and  the  gas  usage  rate  to  increase.  It  is  best  to  have 
smaller  than  the  time  constant  of  the  tracking  miniloop  of  the  mea¬ 
surement  system  by  one  order  of  magnitude.  Generally,  we  should 
let  <  2  seconds. 

2.  The  central  control  circuit  and  the  execution  mechanism 

In  order  for  the  system  to  accomplish  its  tasks  well  in  the 
three  different  working  stages,  as  well  as  to  be  simple,  reliable, 
light  in  weight  and  economic  in  cev;er  consumption,  we  adopted  two 
types  of  switch-controlled  cold-gas  propulsion  mechanisms,  one  with 
a  high  power  and  the  other  a  low  power.  Figure  4  is  the  flow  chart 
of  the  pitch  channel  control  system.  Figure  5  is  the  chart  of  the 
organizational  principle  of  the  execution  mechanism. 

For  convenience,  we  shall  use  the  pitch  channel  as  an  example 
to  describe  the  high  and  low  power  propulsion  systems,  and  we  shall 
consider  the  external  frame  angle  of  the  horizontal  gyro  of  the 
measurement  system  as  approximately  the  pitch  attitude  angle  of 
the  satellite. 
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(1)  The  high  power  propulsion  leading  correction  system:  The 
high  powered  propulsion  leading  correction  system  is  mainly  used  to 
eliminate  the  satellite's  initial  attitude  deviation  stage  and  to 
establish  quickly  the  re-entry  attitude  stage.  In  order  to  increase 
the  resistance  of  the  system,  a  leading  correction  network  is  adopted. 
The  principle  requirement  in  the  design  of  a  high  powered  propulsion 
system  is  to  eliminate  the  initial  attitude  deviation  with  a  minimal 
amount  of  gas  usage  and  to  establish  the  attit ude  needed  for  the 
re-entry  of  the  satellite  with  the  shortest  amount  of  time  and  the 
minimal  amount  of  gas. 

Prom  Figure  4,  it  can  be  seen  that  under  the  condition  that  the 
dynamic  characteristics  of  the  execution  mechanism  and  the  effect  of 
the  time  constant  otK  of  the  leading  correction  network  are  ignored, 
the  law  of  control  of  the  high  powered  propulsion  system  may  be  des¬ 
cribed  with  the  following  equations: 

m=-l,  mA-  -  A, 

0  +  +  (4) 

m=x  +  1  ,mA=  +  A, 

6  +  0'-rK(8+d,)iZ-0D  (  5  ) 

m  =  0.mA^0,  S~(l ~h)90^d 
-^»  +  /C(<)  +  0,)<(l-/,)0o 
(6) 

S  =  mA  =  (  A/,  u-  mFl,)/I  (7; 

Here  0  and  0  are  respectively  program  turn  angle  and  angular 

fc''  - 

speed  in  establishing  the  re-entry  altitude.  A^  and  A  are  respectively 
the  angular  acceleration  produced  by  the  force  moment  of  the  distur¬ 
bance  and  the  control  force  moment  FI , ,  F  is  the  propulsion  force, 

1.  is  the  lever.  Apparently,  during  the  period  of  eliminating  the 

J  • 

initial  altitude  deviation,  0  and  9  are  both  0.  In  oracticai  design, 

P  P 

in  order  to  take  into  consideration  the  effects  of  the  dynamic  char¬ 
acteristics  of  the  execution  mechanism  and  of  the  time  constant  K, 
equations  ( 4  )  — ( 7 )  should  be  modified.  We  shall  briefly  describe  the 
problems  of  selecting  the  principal  parameters  contained  in  the  above 
equations . 
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Figure  5-  Organizational  principle  of  gas  propulsion 
execution  mechanism 

1 — electromagnetic  gate  control;  2 — pitch  jet;  3 — yaw  jet;  3 — roll 
jet;  5 — pressure  adjustment  gate;  6 — gas  bottle;  7 — pressure 
reduction  gate;  3 — 210  atmosphere  0°C;  9 — pressure  adjustment  gate 


Figure  6.  Curves  of  the  transition  time  T  ar.d  total 
propulsion  momentum  1^.  vs.  K. 

1 — transient  oeriod  time  T.  (sec);  2 — total  impulse  I_  (kg/s) 

3— K  (s)  t 

i)  selection  of  the  leading  correction  network  parameter  K: 
Under  the  above  mentioned  simplified  conditions,  the  effect  of  para¬ 
meter  K  on  the  transition  period  time  T^_  and  on  the  total  propulsion 
momentum  may  be  plotted  by  using  either  the  analytical  method  of 
the  simulation  method  as  shown  in  Figure  6.  is  the  time  needed 

to  eliminate  a  certain  amount  of  initial  attitude  deviation.  After 
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the  relative  impulse  of  the  gas  is  given,  the  gas  consumption  rate  Q 
may  be  found  from  the  total  impulse  1^.  by  proportion.  In  establish¬ 
ing  the  re-entry  attitude,  it  corresponds  to  the  time  needed  to  elim¬ 
inate  the  satellite  attitude  angular  deviation  and  the  attitude  angu¬ 
lar  velocity  after  the  pitch  axis  has  turned  to  a  pre-determined  angle 

0  according  to  program.  Figure  6  shows  that  under  the  condition  that 
P 

the  other  parameters  are  relatively  fixed,  and  I  decrease  with  K 
within  a  certain  range,  but  when  K  is  greater  than  four  seconds,  Tt 
even  increases  slightly.  Hence,  it  is  suitable  to  choose  K  =  4  secs. 

ii)  Selection  of  the  switching  characteristic  dead  region  and  the 
stagnant  ring  coefficient:  In  order  to  save  gas,  the  high  powered  propul¬ 
sion  system  will  not  be  used  if  the  attitude  deviation  is  within  ±1°. 

At  the  same  time,  it  is  net  suitable  to  select  too  big  an  0^  if  we 

want  tc  make  the  satellite  pitch  attitude  raoidly  reach  8  +e  (where 

P 

e  is  the  allowed  error).  In  our  design,  we  take  =  1°.  When  the 
stagnant  ring  coefficient  h  is  large,  the  noise  suppression  ability 
is  high  but  the  system  stability  is  not  good.  After  synthetic  consi¬ 
deration,  we  take  h  =  0.1. 

iii)  Selection  of  the  control  of  angular  acceleration  A:  When  the 
effect  of  A  on  T.  and  I-  is  considered,  v/e  can  draw  Figures  7  and  S. 

1  U  I 

From  the  graphs  one  can  see  that  at  K  =  3-6  sec,  when  /}>0.5V,»  ,  both  6 

T.  and  I_  fail  to  decrease  with  increasing  A.  When  A  is  too  small, 

T.,  and  I_  increase  significantly.  Hence,  it  is  suitable  to  select 

,4-0.57... 

iv)  Selection  of  programmed  turn  velocity  6  :  For  simplicity,  reli- 

P 

ability  and  ease  of  implementation,  we  have  chosen  to  use  fixed  pro¬ 
grammed  turn  velocity  0  .  The  selection  of  0  should  be  balanced 

P  P 

between  Tn  (time  needed  to  establish  re-entry  attitude)  and  the  gas 

n 

consumption  rate  Q.  Here  we  try  to  minimize  Q  under  the  condition 
that  T^  ,<  =  60  sec.  In  the  design,  we  chose  0,«»2.67 x  l<T’0//t. 

(2)  Low  powered  propulsion  system:  The  low  powered  propulsion 
system  is  used  chiefly  in  situations  which  requires  highly  precise 
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a^ritude  control  over  long  >  erioric  :f  orbital  motion.  In  order  to 
minimize  the  amount  of  gas  usage  needed  to  satisfy  the  required 
attitude  precision  and  attitude  angular  velocity  with  good  dynamic 
resistance  and  long  life,  we  adopted  a  pseudo-velocity  increment 
feedback  control  system  with  the  discharge  time  constant  Tfd  greater 
than  the  charging  time  constant  T^c ,  as  shown  in  Figure  4. 

The  pseudo- veloc ity  increment  feedback  control  system  is  a  pulse 
regulated  system.  Under  the  situation  of  no  disturbing  force  moment, 
the  system  is  in  an  ideal  limiting  cycle  (1)  operational  state  as 
shown  in  Figure  9(a)  (when  no  speical  measure  is  taken,  it  is  very 
difficult  for  the  system  to  operate  in  the  optimal  limiting  cycle 
(2)  operational  state).  The  equivalent  minimal  propulsion  pulse 
width  T^  output  by  the  system  is 

T-—T  '•>•(>-■ $r) 

where  T  is  the  minimal  culse  width  from  the  pseudo  velocity  incre- 
on 

ment  controller,  AT  is  the  increased  pulse  width  after  the  pulse 
delay  is  taken  into  consideration.  It  is  determined  by  equation  (10) 
where  T  and  T  ,  are  resoectively  the  leading  and  trailing  delay  time 
of  the  electromagnetic  gate.  T^,  and  T  ^  are  respectively  the  lead¬ 
ing  and  trailing  edge  equivalent  time  constants  of  the  propulsion 
jet.  0^  is  the  stagnant  region  of  the  pseudo-velocity  increment  con¬ 
troller.  Here  the  limiting  cycle  angular  velocity  is 


Under  the  condition  that  the  ratio  of  the  satellite  structure  to 
the  jet  gas  volume  is  constant,  we  may  define  the  unit  time  angular 
velocity  increment  to  be  the  weighting  factor  of  the  gas  usage,  denoted 
as  q;  then  during  the  ideal  limiting  cycle  operational  periods  with 
no  disturbing  force  moment, 
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Figure  ? .  The  limit  cycle  operation  of  the 
velocity  increment  feedback  system. 


Here  a  is  the  control  angular  acceleration  produced  by  the 
action  of  the  small  propulsion  f.  T  is  the  period  of  the  limit 
cycle.  It  can  be  seen  from  equations  (11)  and  (12)  that  to  decrease 
9^  and  lower  q,  we  need  to  decrease  ,  and  the  smaller  AT  is,  the 
better.  3ut  in  reality,  AT  is  generally  greater  than  zero,  hence  T, 

u 

is  generally  greater  than  T 

3  on 


However,  under  the  action  of  force  moments,  the  system  is  in  the 

limit  cycle  (1)  and  (2)  operational  state  as  shewn  in  Figure  -o). 
under  most  circumstances ,  it  is  difficult  to  be  in  the  optimal  limit 

cycle  (?)  ccerational  state.  Here  curve  (1)  is  the  phase  trajec¬ 
tory  when  the  disturbing  force  moment  is  less  than  that  for  curve  (2). 
Apparently,  the  smaller  0R  is  designed,  the  easier  it  is  to  cause  pro¬ 
pulsion  on  one  side.  The  bigger  the  disturbing  force  moment,  the 
more  frequent  will  be  the  propulsion.  If  the  Impulse  moment  AHh 
produced  on  the  satellite  by  the  disturbing  force  moment  H^  is  con¬ 
stant,  then  theoretically  speaking,  the  mass  of  the  gas  used  to  over¬ 
come  the  disturbing  force  moment  will  also  be  constant.  If  the  jet 
used  each  time  is  completely  utilised  to  overcome  the  disturbing  force 
moment,  then  the  number  of  propulsions  is 

\M*dt  (/jJ 


From  equation  (13) 
moment  exists,  in  order 
magnetic  gate  operates, 


it  can  be  seen  that  when  a  disturbing  force 
to  reduce  the  number  of  times  the  eleetro- 
to  increase  the  life  of  the  execution  mechanism 
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and  not  to  increase  the  gas  usage,  then  one  should  suitably  increase 
9-  or  T  with  the  condition  that  the  attitude  angular  velocity  for 
fixed  direction  Earth-oriented  observation  should  be  satisfied.  It 
should  be  noticed  that  it  is  not  good  to  increase  Ton  suitably. 

Then  propulsion  from  both  sides  will  again  occur  (see  Figure  9(b) 
curve  4),  causing  the  gas  usage  rate  to  decrease  which  is  not 
desirable  . 


In  summary,  for  the  yaw  channel  which  has  a  large  constant  dis¬ 
turbing  force  moment,  we  have  made  in  our  design  its  limit  cycle 
velocity  about  twice  that  of  the  pitch  and  roll  channels  which  have 
the  smaller  constant  disturbing  force  moment.  The  corresponding  T, 

•  Cl 

satisfies  the  same  kind  of  relationship.  0n  and  T,  are  selected  as 

n  a 

shown  in  Table  1. 


Table  1.  Results  of  eD  and  T,  selection 

„  n  a 


x.  parameter 

channel 

parameter  limit  cycle 
angular  velocity  p,  ! 

equivalent  minimum  propulsion 
pulse  width  T^ 

pitch  and  roll 

yaws 

2.5  x  10_3o/s  | 

5.0  x  10-3o/ 

_ s  1 

50  ms 

100  ms 

Table  2.  Distribution  of  T  and  AT 

on 

-parameter 

channel  — — — __ 

T  (ms) 

on 

A  T  (ms) 

pitch  and  roll 

yaw 

37 

67 

13 

13 

It  may  be  seen  from  equation  (11)  that  9„  is  determined  by  T . . 

•  n  Cl 

In  selecting  9^  and  T^,  we  should  consider  the  problem  of  selecting 
a.  In  the  design  we  have  chosen  a  =  0.1  °/s^.  It  is  much  greater 
than  the  angular  velocity  produced  by  the  disturbing  force  moment, 
hence  making  sure  that  the  high  and  low  powered  propulsions  will  oper¬ 
ate  in  conjunction.  At  the  same  time,  a  is  not  too  small  so  that  it 
is  easily  implemented. 
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In  addition,  we  may  suitably  distribute  Ton  and  AT  from 
equation  (8)  in  accordance  with  necessary  analysis  and  experimenta¬ 
tion.  In  our  design,  the  distribution  of  TQn  and  AT  is  shown  in 
Figure  2. 

As  is  well  known,  as  far  as  pseudo-velocity  increment  feedback 

control  system  is  concerned,  when  the  attitude  deviation  9  reaches  a 

certain  value  6  ,  the  system  will  be  in  a  continuous  propulsion  state, 
s 

We  call  the  attitude  deviation  9  which  starts  such  a  state  by  the 
name  "saturated  jet  angle".  It  may  be  expressed  as 

9<  =  Kia+{\-h)6d  8 

It  is  an  important  parameter  to  guarantee  that  the  high  and  low 
powered  propulsions  operate  in  conjunction.  In  our  design,  we  have 
considered  the  fact  that  the  high  powered  propulsion  switch  charac¬ 
teristic  stagnant  region  9^  is  1°.  To  guarantee  that  the  high  and 
low  powered  propulsions  do  operate  in  conjunction,  we  took  6s  to  be 
0.6-0.75°. 


Besides,  under  the  condition  that  the  operational  characteristics 
of  the  execution  mechanism  and  that  of  the  measurement  system  are  not 
considered,  the  system  may  be  roughly  seen  as  a  second  order  linear 
system  with  resistance  ratio  ;  and  natural  frequency  by  the  des¬ 
criptive  function  method  . 


On 
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After  the  basic  relationship  and  data  are  given  as  above,  we 
can  then  proceed  to  make  selections  on  some  of  the  other  major  para¬ 
meters  of  the  low  powered  propulsion  system  without  much  difficulty. 


i)  Selection  of  the  switching  characteristics:  The  stagnant 
region  9^  should  match  the  measurement  error  of  the  measurement 
system  in  order  to  satisfy  the  attitude  precision  requirement  of  the 
orbital  motion  while  the  switch  stagnant  region  of  the  low  powered 
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propulsion  system  8^  should  be  much  smaller  than  that  of  the  high 
powered  propulsion  system  9^  so  that  the  low  powered  system  is  the 
one  that  operates  in  the  orbital  motion.  We  chose  8^  =  0.3°. 

Stagnant  cycle  coefficient  h  affects  TQn  on  one  hand, while  it 
also  affects  to  a  certain  extent  the  ability  of  the  system  to  suppress 
noise.  From  equation  (9),  we  know  that  when  the  other  parameters  are 
held  constant,  the  smaller  h  is,  the  smaller  TQn  will  be,  while  the 
ability  to  suppress  noise  will  be  weaker.  Taking  into  consideration 
that  the  selection  of  the  other  parameters  has  already  provided  the 
system  with  a  strong  noise  suppression  capability,  we  chose  h  =  0.1 
after  some  simulation  testing  and  analysis. 

ii)  Selection  of  the  feedback  network  parameter:  In  order  to 
make  sure  that  the  high  and  low  powered  propulsion  work  in  conjunction, 
we  take  8  =  0.7.  From  equation  (14),  it  may  be  determined  that  Kp  = 

4 .  3  secs  . 

In  order  to  make  sure  that  the  pulse  width  Tqr  is  minimal,  we 
may  compute  from  equation  (9)  the  charging  time  constant  Tfc  based 
on  the  other  parameters  already  determined.  In  order  to  make  sure 
that  the  system  possesses  a  large  resistance,  we  can  determine  the 
discharging  time  constant  from  equation  (15).  In  the  design,  we 
take  ;  =  1 .  Thus,  the  results  of  selecting  Tfc  and  Tf^  are  shown  in 
Table  3- 


The  results  of  the  selection  of  the  principal  oarameters  for  the 
high  and  low  powered  propulsion  systems  may  easily  be  transformed 
into  the  requirements  on  the  dynamic  characteristics  and  the  sice  of 
power  of  the  propulsion  execution  mechanism.  The  normalised  values 
of  the  high  and  low  propulsion  strengths  of  the  various  channels  are 
shown  in  Table  4. 


In  the  preliminary  system  design  stage,  a  large  amount  of  math¬ 
ematical  simulation  experimentation  was  carried  out.  After  the  ini¬ 
tial  production  samples  were  produced,  experiments  with  the  real  and 
semi-real  objects  were  done  to  verify  the  system  design.  The  j  oi.nt  ocerat 
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Table  3.  Results  on  the  selection  of  T-  and  T 


Table  4.  Normalised  values  of  the  propulsion  strength 
for  all  channels 


— type 

channel  ■ — ^ 

high  power 
(gm) 

low  power 
(gm) 

pitch 

1000 

200 

yaw 

1000 

200 

roll  i 

350 

70 

of  the  high  and  low  powered  propulsion  in  each  channel  and  the  effect 
on  the  coupling  of  the  three  axes  are  also  studied.  Based  on  that 
study,  the  prototype  was  developed.  The  interface  relationship 
between  the  main  system  and  its  various  sub-systems  and  the  require¬ 
ments  of  them  were  then  finalised.  In  ‘•he  development  cf  the  proto¬ 


type,  a  large  amount  :f 

satisfy  the  reliability 
shall  not  describe  them 


an  a  ly  s  i  s  an  d  e  reel  me  r.  r 

required  for  the  satell 
owing  to  the  limitation 


it  ion  was,  carried  cut 
ce  after  launch.  Vie 

of  space. 


III.  RESULTS  OF  FLIGHT  TESTS 


In  order  to  examine  the  operational  characteristics  of  the  atti¬ 
tude  control  system  in  the  flight  test  and  to  obtain  some  useful 
data,  stellar  cameras  that  are  capable  of  discriminating  the  preci¬ 
sion  of  the  attitude  control  system  are  installed  in  the  satellite 
besides  setting  up  some  remote  sensing  parameters  related  to  the 
attitude  control  system.  We  shall  briefly  describe  below  the  oper¬ 
ational  status  and  preliminary  analytical  results  of  the  attitude 
control  system  in  three  different  operational  stages  according  to  the 
data  obtained  in  the  three  flight  tests. 
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1.  The  stage  during  which  the  initial  attitude  deviation  is 
being  eliminated. 


Remote  data  indicate  that  among  the  three  satellites,  the  first 
one  launched  had  a  relative  bad  initial  attitude  deviation.  The  two 
satellites  launched  later  had  made  great  improvement.  With  the  first 
satellite  as  an  example,  Table  5  shows  its  initial  attitude  devia¬ 
tions.  The  dynamic  procedure  to  eliminate  the  initial  attitude 
deviation  is  shown  in  Figure  10. 


Table  5.  The  initial  attitude  deviation  of  the  first 
satellite 


attitude 

|  initial  attitude 

initial  attitude 

— deviation 

cnannel 

angle(degrte ) 

angular  velocity 

(Vs) 

pitch 

+1.13° 

-2.24 

yaw 

+1.63° 

+0.17 

roll 

+0.98° 

-0.17 

2 

¥  ' 
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Figure  10.  Dynamic  process  in  the  elimination  of  the 
initial  attitude  deviation  for  the  first  satellite, 
(from  remote  gyro  data). 

1 — roll  angle;  2 — yaw  angle;  3 — pitch  angle;  4 — time  (sec) 
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The  operational  properties  of  the  attitude  control  system  may 
be  seen  from  Table  6. 


TABLE  6 

"---^content 

satellite  n3'>--. 

max.  angle  of 
adjustment  or  procedure 
[  characteristic  point 

time  on  entering 
0.7° 

i 

-5-36° 

12  sec 

II 

single  adjustment  decay 

12  sec 

III 

single  adjustment  decay 

|  6  sec 

The  above  situation  indicates  that  the  operation  of  the  atti¬ 
tude  control  system  during  the  initial  attitude  deviation  elimination 
stage  is  normal  and  is  very  close  to  the  theoretical  analysis.  The 
system  has  a  strong  initial  attitude  capture  capability. 

i  2.  During  the  orbiting  stage 

(1)  Concerning  the  accuracy  of  the  attitude  control.  The  actual 
accuracy  attained  by  the  attitude  control  system  is  obtained  accord¬ 
ing  to  the  stellar  camera  data.  The  stellar  camera  obtains  first  the 
attitude  angle  of  the  satellite  relative  to  the  Earth-centered  astro¬ 
nomical  coordinate  system  by  using  the  method  of  photographing  the 
stellar  sky,  and  then  transforms  it  into  the  attitude  angle  relative 
to  the  Earth-centered  orbital  coordinate  system  from  the  orbit  track¬ 
ing  data.  The  satellite  attitude  obtained  from  the  stellar  camera 
data  is  more  accurate.  Figure  11  shows  the  typical  attitude  data 
obtained  from  the  stellar  camera  data.  The  attitude  angle  and  atti¬ 
tude  angular  velocity  deviation  actually  attained  are  displayed  in 
Figure  7- 


TABLE  7.  Attitude  deviation  actually  attained 


— ^2_ontent 

attitude  angle 

attitude  angular 

channel^'"' — ^ 

deviation 

deviation 

pitch 

/a 

II 

o 

—0 

«  =  0.01 

roll 

.<  =  0.7 

II 

O 

O 

ro 

yaw 

^  =  1.5 

v<  =  0.02 

11 


— 

*» 
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if 

«*<»> 

Figure  11.  Typical  attitude  data  from  stellar  camera 
1 — yaw  angle;  2 — pitch  angle;  3 — roll  angle;  ii  —  time  (sec) 

(2)  Concerning  the  spectral  density  of  the  infrared  levelmeter 
noise:  Preliminary  analysis  has  been  carried  out  on  the  infrared 

levelmeter  output  noise  based  on  the  remote  data  on  the  infrared 
levelmeter  and  the  gyro  frame  angle  during  the  orlbtal  motion  stage. 
In  the  analysis,  a  smooth  random  model  is  assumed  for  the  noise.  The 
power  spectral  density  of  the  infrared  levelmeter  output  noise  then 
can  be  calculated  through  Fourier  transform.  We  used  a  method  simi¬ 
lar  to  that  in  [U],  After  making  a  dynamic  correction  on  the  measur¬ 
ing  system,  it  is  found  that  the  band  width  of  infrared  levelmeter 
output  noise  power  spectral  density  is  approximately  12  times  the 
orbital  angular  velocity.  Because  there  is  not  enough  continuous, 
analysabie  data,  this  is  only  a  preliminary  result  usable  for 
reference . 


Figure  12.  Typical  operational  character  in  establishing 
re-entry  attitude  (obtained  from  gyro  remote  data'' 

1 — output  from  horizontal  gyro  large  angle  transducer;  2  —  comnand 
time;  3 — output  from  horizontal  gyro  small  angle  transducer 

3.  During  the  stage  of  establishing  re-entry  attitude 

During  the  re-entry  attitude  establishment  stage,  the  attitude 
control  basically  operates  in  the  same  way  in  the  three  satellites. 
Its  typical  operational  character  is  shown  in  Figure  12.  It  shows 
the  complete  dynamic  process  in  establishing  re-entry  attitude  by 
the  attitude  control  system  of  the  third  satellite. 

From  Figure  12,  it  may  be  seen  that 

(1)  the  output  from  the  horizontal  gyro's  outer  frame  large 
angle  transducer  indicates  that  the  satellite  had  rotated  about  the 
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pitch  axis  in  the  negative  direction  through  an  angle  0^  degrees. 

Its  dynamic  process  agrees  exactly  with  the  theoretical  analysis. 

(2)  The  output  from  the  horizontal  gyro's  outer  frame  small 
angle  transducer  indicates  that  at  the  beginning  of  attitude  adjust¬ 
ment,  the  output  rises  sharply  with  a  positive  sign  to  reach  a  maxi¬ 
mum  over-adjustment  angle.  After  a  certain  time  period  T  ,  the  output 
is  almost  zero.  It  can  be  ascertained  from  this  process  chat  the 
program  mechanism  turns  through  0  degrees  with  speed  6  according 
to  a  pre-determined  program  in  a  pre-determined  direction.  The 
satellite  will  also  attain  this  angular  velocity.  Finally,  the 
satellite  is  stabilized  in  the  range  9^  +  e .  The  complete  attitude 
adjustment  process  takes  about  55  seconds. 

The  above  indicates  that  the  process  of  establishing  re-entry 
attitude  basically  agrees  with  the  prediction.  The  operation  is 
entirely  normal. 

•4.  Concerning  the  gas  usage  rate  of  the  attitude 
control  system 

The  experimental  results  of  three  flights  indicate  that  * r-  gas 
usage  rate  of  the  attitude  control  system  (whether  theoretics  esti¬ 
mate  or  actual  consumption)  gradually  decreased,  especially  for 
the  third  satellite.  Its  gas  consumption  is  two  times  smaller  than 
that  of  the  first  satellite.  The  main  reasons  for  this  decrease  are: 

(1)  based  on  the  experiment  of  the  first  two  flights,  we  have 
modified  the  parameters  of  the  attitude  control  system  of  the  third 
satellite,  especially  the  principal  parameters  of  the  measuring  system. 
We  have  also  improved  the  noise  suppression  ability  of  the  system  and 
reduced  gas  consumption  rate. 

(2)  we  have  added  a  front  filter  to  the  infrared  levelmeter 
optical  system,  improved  the  bandwidth  of  the  meter  and  reduced  the 
effect  of  the  sun,  the  moon  and  the  cold  cloud  on  the  output  of  the 
infrared  levelmeter,  hence  reduced  the  gas  consumption. 
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CONCLUSIONS 


xv 


The  results  of  the  three  flight  tests  indicate  that  the  planning 
and  design  of  the  satellite  attitude  control  system  of  our  first 
generation  re-entry  type  Earth-oriented  fixed  direction  observation 
technique  is  correct.  The  system  worked  reliably  and  operated 
normally.  In  particular,  the  improvement  on  the  third  satellite  is 
effective.  The  mission  is  completed  satisfactori ly  and  some  useful 
data  have  been  obtained.  Some  experience  has  been  accumulated  which 
has  established  a  foundation  for  developing  similar  satellite  control 
systems  in  the  future. 

In  the  future,  it  is  necessary  to  improve  on  the  accuracy  of 
the  attitude  control  system,  especially  the  yaw  attitude  accuracy. 

The  operational  lifetime  of  the  satellite  should  be  lengthened  with 
reduced  development  budget.  System  engineering  and  modern  control 
theory  should  be  applied  well  to  the  design  and  development  of  atti¬ 
tude  control  systems. 
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GUIDANCE  ACCURACY  DETERMINATION  BASED  ON  DATA  MEASURED 
EXTERNALLY  DURING  PROPULSION-FREE  STAGE 


Zhu  Wenxuan 

ABSTRACT 

We  propose  in  this  paper  to  measure  the  position 
vector  of  a  space  vehicle  in  its  propulsion-free  orbit 
with  radar  and  compare  it  with  the  normalized  oribtal 
position  vector.  By  treating  the  difference  variationally 
with  matrix  theory,  it  is  possible  to  determine  the  total 
guidance  aberration  of  the  guidance  system  instruments  on 
the  space  vehicle  when  the  thrust  ends.  This  may  provide 
an  auxiliary  method  for  accuracy  checking  to  compensate 
for  the  shortcomings  of  the  ballistic  measurement  devices 
during  propulsion  as  well  as  the  accurate  position  vector 
of  the  normalized  orbit  during  atmosphere  re-entry. 

I .  INTRODUCTION 

The  accuracy  of  the  guidance  system  and  inertial  system  of  car¬ 
rier  rocket-type  space  vehicles,  especially  the  accuracy  In  the 
measurement  of  inertial  instruments,  is  a  problem  of  common  concern 
for  departments  designing,  testing  and  using  space  vehicles.  One  cf 
the  important  goals  of  developing  highly  accurate  remote  measurement 
systems  and  large,  complex  and  highly  accurate  ballistic  measurement 
systems  Is  to  determine  the  accuracy  of  the  guidance  system  and  its 
inertial  sensitivity  measurement  instruments,  as  well  as  to  determine 
the  reliability  of  the  ground  test  accuracy  of  the  inertial  instru¬ 
ments  . 

When  the  space  vehicle  is  under  propulsion  (engine  thrust),  there 
exist  large  acceleration  and  sudden  impulse  ,  random  and  large  rota¬ 
tions  about  the  center,  complex  and  severe  vibrations  of  the  three 
axes,  large  variations  In  pressure  and  temperature  as  well  as  all 
kinds  of  radiation  interference.  For  inertial  instruments  working 
in  this  kind  of  environment  under  conditions  far  worse  than  those  in 


the  laboratory  or  ground  level  factory,  it  is  necessary  to  use 
vehicle  ballistic  measurement  parameters  one  order  of  magnitude  more 
accurate  than  those  measured  by  the  inertial  instruments  to  compare 
with  the  internal  remote  measurement  parameters,  and  to  treat  and 
analyze  the  huge  amount  of  accurate,  reliable  data  in  order  to 
obtain  reliable  results  in  determining  the  operational  accuracy,  in 
measuring  the  sensitivity  and  in  varying  the  parameters  being  mea¬ 
sured.  However,  for  certain  types  of  vehicles,  it  is  often  impossibl 
to  realize  the  required  accuracy  due  to  limitations  in  material  condi 
tions  and  measurement  accuracy. 

In  the  known  gravitational  field  of  the  Earth,  satellite  type 
vehicles  move  in  thrustless  free  flight  state.  By  tracking  these 
vehicles  with  ballistic  measurement  systems  and  measuring  their 
orbital  parameters  as  well  as  analyzing  and  carrying  out  computation¬ 
al  treatments  on  the  measured  data,  it  is  possible  to  obtain  a  synthe 
tic  determination ,' including  both  propelled  and  propulsion  free 
flights)  for  the  total  guidance  accuracy  of  the  inertial  system,  and 
to  compensate  for  the  shortcomings  in  determining  the  accuracy  of 
the  inertial  system  In  the  propelled  stage.  In  this  paper  we  shall 
discuss  the  theoretical  aspect  of  this  method. 

II.  THE  BASIC  PRINCIPLES 

To  determine  the  orbit  of  a  space  vehicle  with  the  externally 
measured  parameters  of  the  propulsion- free  (0  thrust)  vehicle  is  a 
common  method  used  by  departments  designing,  testing  and  utilizing 
the  vehicle.  Using  these  externally  measured  parameters  to  determine 
the  total  guidance  aberration  (accuracy)  is  a  new  direction  in  using 
these  parameters. 

After  the  thrust  on  the  vehicle  has  ended  (engine  shut-off), 
the  velocity  of  the  vehicle  may  reach  thousands  of  meters  per  second 
at  a  height  of  more  than  100  kilometers.  Over  most  of  the  flight 
time,  the  atmosphere  in  the  orbit  is  very  rare  so  that  aerodynamics 
may  be  neglected.  In  comparison  with  the  gravity  of  the  Earth,  one 
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nay  also  neglect  the  effect  of  photon  pressure,  radiation  pressure 
and  gravitational  attraction  of  other  heavenly  bodies  on  the  vehicle. 
Hence,  the  vehicle  travels  along  a  pseudoelliptic  orbit  solely  under 
the  Earth's  gravitational  attraction  (if  we  regard  the  Earth  as  a 
sphere  and  not  a  spheroid,  i.e.,  we  neglect  the  oblateness  of  the 
Earth,  then  the  gravitational  force  is  an  inverse  square,  central 
force.  The  venlcle  will  move  in  a  plane  along  a  pseudoelliptical 
orbit).  The  principal  parameters  of  this  pseudceliipse  such  as  the 
semi-ma.jcr  axis  and  eccentricity  are  uniquely  determined  by  rhe  velo¬ 
city  i  f(V..Vr  .V,)  and  position  rf(.*.y.s)  when  the  thrust  is  ended.  The 
pseuaoellipse  of  satellites  ana  spacecraft  does  net  intersect  with 
the  Earth.  They  will  rotate  around  the  Earth.  The  pseudoellipse  of 
a  guided  missile  and  its  warhead  does  intersect  with  the  Earth.  It 
has  to  re-enter  the  atmosphere  and  fall  in  a  pre-determined  target 
region.  The  positions  of  at  least  three  points,  ,  on  the 

orbit  of  the  vehicle  in  a  propulsion- free  region,  should  be  measured 
if  one  uses  external  ballistic  position  measurement  systems,  such  as 
single  pulsed  radar  or  photographic  theodolites  with  lower  accuracy 
than  that  used  in  propelled  flight .  The  angle  subtended  by  each 
pair  of  points  should  be  greater  than  G.5  ^  1°.  The  true  pseudo- 
elliptic  orbit  of  the  vehicle  may  then  be  obtained  by  applying  proper 
computations  on  Y3>  This  is  then  the  ordinary  orbit  determination 
method. 


The  velocity  pi  and  position  fr.  at  the  thrust  termination 
point  (engine  shut-off  point)  may  also  be  derived  through  proper  com¬ 
putations  from  the  y3  on  two  points  which  subtend  an  angle  at  the 
center  of  the  Earth  greater  than  0.5  ^  1°.  The  combination  of  P'.r’ 
is  needed  to  satisfy  the  guidance  condition  or  for  the  guidance  system 
to  fulfill  its  mission.  Their  deviation  from  the  theoretically  com- 

—  r* 

puted  velocity  p’t  and  position  at  ?(  at  tA  when  the  thrust  stops: 

iP’.n-r:  (2.!, 

dr  I  **  j*—  r  t 

are  the  basic  data  for  calculating  the  total  guidance  error.  The 
deviation  from  the  target  of  the  missile  impact  point  calculated  from 
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1 — propelled  stage;  2 — ascending  stage;  3 — standard  orbit; 
i — descending  stage;  5 — re-entry  stage;  6 — center  of  Earth 

dpi  and  dr'  (off-target  amount)  and  the  deviation  of  the  number  25 
oribts  of  the  satellite  or  spacecraft  all  belong  to  the  total  guid¬ 
ance  error.  After  reducing  from  if,,  if?  ,  the  errors  dp'.Cir'  due 
to  guidance  principle  (method),  the  remainder  is  then  the  sum  of  the 
error  due  to  the  guidance  system  instruments. 

For  re-entry  missiles,  the  normalized  ballistic?*.-  for  the 
atmosphere  re-entry  stage  may  be  established  from  1  ''<?'.  deter¬ 
mined  from  the  externally  measured  data  ?•  from  the  real  orbit.  If 
an  external  ballistic  tracking  measurement  is  carried  out  on  the 
position  of  the  re-entry  stage  of  the  missile,  then  the  deviation 

ar. (2-2) 
is  the  so-called  re-entry  error  due  to  the  aerodynamic  disturbance 
of  re-entry.  Or  the  re-entry  error  is  the  difference  between  the 
deviation  of  the  impact  point  <5 Qr  calculated  from  Ji  \  '  and 
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the  deviation  bQ, 


cf  the  actual  point  of  impact. 


We  shall  demonstrate  the  above  principles  with  a  re-entry 
missile  as  an  example.  Its  orbit  is  shown  in  the  above  diagram. 

III.  DYNAMIC  EQUATION  FOR  THE  PROPULSION-FREE  STAGE 

We  shall  establish  our  dynamic  equations  in  the  vertical  coor¬ 
dinate  system  Oxyz  with  the  launching  point  0  as  the  origin  with  oxz 
in  the  horizontal  plane  through  0  and  ox  basically  along  the  flight 
direction.  Oxyz  rotates  with  the  Earth.  The  propulsion-free  stage 
(free  flight  stage)  starts  from  the  thrust  ending  at  point  K, 
through  the  ascending  stage,  the  descending  stage,  at  the  re-entry 
stage  and  ends  with  the  point  1  on  Earth.  The  vehicle  is  acted  on 
by  gravitational  force,  centrifugal  force,  Copernican  force  and 
aerodynamic  force.  The  propulsion- free  dynamic  equations  are 
obtained  after  the  aerodynamic  lift  force  and  side  force: 

initial  condition  tf  =  tb 

o 


y 

i-v. 


,>  =  JL  r  p'Sf  ( p 


2  J  V, 


(3-D 


where  ,  3^,  Cf  are  respectively  the  mass,  reference  cross-sec¬ 
tional  area  and  aerodynamic  frictional  coefficient  during  the  pro¬ 
pulsion-free  stage,  x,  y,  z  and  Vv,  V  .  V ^  are  the  position  and 

x  y  ^ 

velocity  components  of  the  vehicle  in  Oxyz. 

The  linking  acceleration  components  are: 
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-(*  +  *,)];  j  ^3-2 

where  to  is  the  rotational  angular  velocity  of  Earth  and  the  flight 
direction  coefficients  are 


b,  *«  cos  W ,  cos  Am 
b,~  sin  W, 
b.*  —  cos  W,  sin  Am 


(3-3) 


W.  is  the  astronomical  longitude  of  the  point  0.  A  is  the  astrono- 
t  m 

mical  position  angle  of  ox.  The  components  of  the  radial  distance 
EO  of  the  point  0  to  the  center  of  the  Earth  in  Oxyz  are: 


x,  —  /?,.  sin  )it  cos  Am  | 

y,  =  A\,cosM,  > 

sin  y.sin  Am  1 


(3-4) 


H  is  the  height  of  point  0.  The  radial  distance  of  any  point  on 
the  Earth's  elliDsoid  of  rotation  is 


Earth  center  longitude: 


/?=  i  — 

v  1  —  e'cos  :o. 


.  .  _,c 

^,=  sin  '  — 


The  equator  surface  height: 

£=*bjt  +  b,\'  +  b,s  +  R„  sin 


(3-5) 


(3-6) 


(3-7) 


The  radial  distance  of  the  vehicle  position: 

r-V(x  +  x,)'+(v- v.J't  (*  +  *,)*  (3-8) 

The  difference  y  between  the  geographical  longitude  and  the  Earth 
center  longitude  of  a  point  on  Earth  surface  is  calculated  by: 

sin /s —a*  tin  2$, =s=  a*  sin  2^.  (3-9) 
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The  geographical  longitude  at  the  launch  point  is  taken  as  W  =  . 


Equation  (3-9)  is  used  to  compute  p  , 
The  relative  velocity  of  the  venicle 


V  =  JV\  +  Vl  +  l-'. 

The  coperican  acceleration  components  are 

V',~2a(bJ\-bJ',).  j> 
t..=2o>(by,-b,r,).  j 


( 3-10 ) 
( 3-ID 

(3-12) 


The  gravitational  acceleration  of  an  ellipsoid: 

R!  .  ,R.  -  ■** 


»,  k;  .  k.  /  .  ±  \  ! 

+JT‘\  3  r J  ~  1 )_) : 

(v-*#  y,  *;  *,  =  *  +  *,) 


where  is  the  equatorial  radius  of  Earth,  gmR\—GM 


(3-13) 

is  the  gravi¬ 


tational  constant  of  Earth. 


ad  is  the  oblateness  of  Earth. 


With  the  given  initial  conditions  at  ,  the  dynamic 

equations  (3-1)  are  solved  to  obtain  the  diagram  in  the  last  section, 
called  the  normalized  or  reference  orbit.  From  the  normalized  orbit, 
the  velocity  ?(*'..  V,.  V .)  and  position  Hx.  y.  z).  of  the  vehicle  at 
any  arbitrary  time  may  be  obtained. 

IV.  STATE  TRANSFORMATION  MATRIX 


At  t'—ti  ,  the  orbital  parameter  deviations  are  given.  For 
example , 

AK>l<r",  AK'-10-'\  AP'-KT". 

A*'- 1000"  Ay- 1000"  Az'=1000"; 

Taking  one  of  these  at  a  time,  we  can  calculate  6  deviation  orbits 
and  also  the  state  transformation  matrix  at  each  time  instant  on  the 
orbit  =  U .  For  example,  at  the  point  1  (i  =  l)  at  t  =  t ^ : 

(4-1) 

LV,  6.',  bi.  b,'t  6,  J 

At  point  (t  =  t  ),  it  is 
m 


For 

example , 

at 

4- 

O 

he  po 

K  *>!> 

bi. 

bi. 

£>; 

b.!,  b‘ ,  bi, 

bi. 

bi , 

K  | 

b}>  bi.  bJ, 

bi. 

bi. 

6.;  J 
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f 


~K  b- 

b-  b‘ 

far,  far. 

b~  6* 

b: ;  fa- 

b-  b-  j 

.b‘  6  r. 

6;,  far. 

/■  »  /s"  1 

(4-2) 

d* 

jr*  —  jr  ‘ 

Ax 

*  ~  d/J  ~ 

AA/ 

-'a  ?J  * 

ay 

k~  dVm 

xirPL 

a;.' 

II 

(4-3) 

d2 

■»  •  —  z  • 

Ai* 

k~dkt~ 

aa/ 

A/T  ; 

A- 1.2, 3. 4. 5. 6, 

A'-  (A*).,,  =  [A,  A,  A,  A,  A,  A.]r  -  [K  (  VI  V '  x'  y /  giy , 


x‘ >  x‘  are  the  x  values  of  point  1  at  t=t.  on  the  normalised  orbit 

and  the  deviation  orbit;  correspondingly,  y> ,  y-  and  have 

the  same  meanings.  Hence,  the  orbital  parameter  deviation  vector 
or  column  matrix  expressed  with  the  transformation  matrix  is: 

A/?'«5,AA/  (4-4) 


where 


Ar Ar|]T  »  [Ax1  Ay  Ar']r . 


For  example,  at  point  i  *  1 


rArj  i 

rAx' 

I 

; 

A/?' 5= 

Ar'  L 

1 

Av' 

'  1 

J  — 

.  J 

-  A*»  J 

. 

fa/. 

fa:, 


b,1, 

fa,',  fa/. 

fa/. 

•  Ar; 

fa.', 

6:’.  fa  A 

fa/. 

!  af; 

fa,', 

fa  .  ,  fa 

fa’. 

|:  Ax/ 

Av'  i 

I 

k  Ai'  j 

The  same  form  holds  at  i  *  m.  We  only  need  to  chance 
script  1  into  m.  At  points  1  and  m,  we  have: 


(4-5) 


(4-6) 


the  super- 


^  Ar  [  n 

(  Ax'  N 

( h ■' 

bi, 

fa/. 

fa/. 

fa/, 

fa/. 

,ak  : 

i  Ar1,  j 

Ay' 

!  fa/. 

b-!. 

fa/, 

fa/. 

fa/, 

fa,'. 

i  j  AK  J 

Ar  (  i 

Az' 

!  !  fa,', 

fa,', 

fa/. 

fa,'. 

fa/, 

fa/. 

Ar;  . 

Ar?  j 

1  Ax" 

!  ;fa- 

fa,": 

far, 

far. 

far. 

far. 

ax/  i 

i  ArT  | 

Ay" 

1  l6’“ 

fas 

far, 

far. 

far. 

fa; 

l!  Ay/  : 

k  ArT  J 

kA zm  / 

far, 

far, 

fa: 

far. 

fa; 

-1  k  A 2<  / 
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This  may  be  expressed  with  sub-matrices  as  follows 

_  f  Ai?'  ~| 

A/?'*  =  „  1 

LAi?*- 


(4-8) 


-Bm- 


(4-9) 


A f  as  shown  in  (4-5).  as  shown  in  (4-1)  or  (4-2). 
The  sub-matrices  3,  and  3  are  both  3x6  rectangular  matrices.  The  sub- 
matrices  ARX  and  AR'71  are  born  3x1  row  matrices.  Hence 


\R-m  =  B,mW 

f 

where  AX*  is  still  a  6x1  row  matrix. 


A,}/-  [A-i,  Ax,  Ax,  Ax,  _U,  A/Jr  =  [A^(  AFi  \l".  Ax' Ay/ A  s'Y 


(4-10) 


(4-11) 


The  super-script  T  outside  the  square  bracket  means  the  transpose  of 
the  matrix. 


V.  CALCULATION  OF  THE  ORBITAL  PARAMETER  DEVIATIONS 
AT  THE  THRUST  ENDING  TIME  tx 


The  position  vector  components  of  the  vehicle  at  t  =  t^  during 
the  propulsion-free  stage  are  measured  with  radar  or  an  optical 
system.  The  vector  or  matrix  is 

$•  =  /?;  +  «'  (5-1) 


where 


R\=*lx\  y\  r', ]T 


(5-2) 


is  the  actual  position  of  the  vehicle  at  point  i  during  the  propul¬ 
sion-free  stage,  v1  is  the  random  measurement  error  of  the  measuring 


instrument 


r'rV.y 


(5-3) 

Generally,  the  average  or  mathematical  expectation  of  the  measure¬ 
ment  error  is  assumed  to  be  zero.  The  difference  between  the  compo¬ 
nents  of  the  actual  orbit  and  the  normalized  orbit  is 


Ax;  -x;— ) 

A y',-y',-y  j> 


(5-4) 


simplified  in  matrix  form  as 


Atf',  =*(Ax;  Av;  Ax’,  ]? 


(5-5) 
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;o  ( 4—6 ) 


A/v  =B, ,AA' 


(5-6) 


The  deviation  of  the  measured  value  l'~  is 

AA‘  =  A/e-,4.v'  =  fi.AA/+v‘  (5-7) 

f 

Now  AX  is  the  unknown  to  be  found.  Expressed  by  AX^,  we  get 

(5-3) 

AX,  =[AF.,  At',,  AK„  Ax,  Ay,  Ar,]T  (5-9) 


If  there  is  only  one  observation  point  1,  then 

fl,AX,=*A  S'-v' 


(5-10) 


It  is  an  equation  system  with  3  equations  but  6  knowns .  AX^  is 
unaetermin  d . 


If  there  are  2  points  of  observation  1,  m,  i.e.,  observed  values 

at  two  times  t.  and  t  ,  then 

1  m 

5, .AX,  =»AS'*-»'*  =  A/?'*  (5-11) 

Its  equation  number  and  unknown  number  are  both  6.  When  the  angle 
of  the  two  points  i,  m  subtended  at  the  center  of  the  Earth  is  large 
enough,  i.e.,  larger  than  1  degree,  the  matrix  is  non-singular. 

Hence,  the  AX^  to  be  found  is  definitely 

AA, (5-12) 
It  is  the  deviation  o^  the  orbit  parameter  that  satisfies  the  guid¬ 
ance  condition  at  the  moment  when  the  thrust  ends  and  that  we  desire 
to  calculate  in  the  section  "Discussion  of  the  Basic  Principles". 
Velocity  deviation  Is 

Ap/-r/—  r'-IAt'.' At''  At'']’  =CAk,l\kriAt',l]  (5-13) 

and  position  deviation  is 

Ar'  =  r [Ax-'  Ay'  A*']'  =■  [  \x,'Ay  A.’,]'  (5-lA) 

When  v'%o,  v"*o  ,  the  value  found  from  (5-8)  contains  errors,  and 
is  not  the  true  value  of  AX1 . 


If  the  number  of  observation  points  N  Is  larger  ,..an  2,  then  the 
data  for  any  two  distinct  points  1,  m  may  be  used  to  calculate  a  set 

J  f 

of  AX“  with  (5-8).  Then  a  better  estimated  value  for  AX^  is 


AL  =*-rS 

V  t  *  1 


(5-15) 
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where  p  >  2  is  the  number  of  observation  points. 


If  there  is  a  large  number  of  .'cserv'ti  or.  points  S,  i  =  1,  2... 
k .  .  .  1 .  .  .m.  .  .  r. .  .  .  N,  then  after  numerical  treatment,  we  can  get  the 

r* 

best  estimated  value  of  AX1  to  be  AX~.  For  each  i  point,  equation 

i 

(5-7)  holds.  The  AS  of  N  points  has  the  same  form 

AS  =  5AA.  +  v  (5-16) 

where 

AS~[aS'  AS-*— as*—  AS'  — AS'.-AS’—AS"]’  (5-17) 

is  a  3N  x  1  row  matrix.  The  submatrices  AS""  are  3x1  row  matrices. 
The  transformation  matrix  B  is  formed  from  3x6  rectangular  matrices 


(5-18) 


Apparently,  b  is  a  31'  x  6  matrix.  3  are  still  computed  from 
(4-1)  or  (  A  -  2  )  and  (4-3).  AA.-CAP.,  AK,,  A*',,  Ax;  A„v,  A*',]t 
is  a  6  x  1  row  matrix.  When  N  >  2,  equation  (5-16)  must  be  over¬ 
determined  after  expansion,  i.e.,  the  number  of  equations  3N  is  far 
greater  than  6,  the  number  of  unknowns  AX^.  The  order  of  the  matrix 
3  =  row  number  *  6.  (The  so-called  order  is  the  large  order  of  the 
non-zero  sub-matrix  of  B).  According  to  matrix  theory,  it  can  be 
proved  that 

BTB=F  (5-19) 

is  a  non-singular  symmetric  matrix  of  order  6,  i.e.,  the  determinant 

\F\-\BTB\-Q  (5-20) 

Hence,  the  inverse  matrix 

F-'-WBT'  ( 5-21 ) 

T 

exists,  and  is  the  only  matrix  of  order  6.  Obviously,  E  is  a  6  x 
3N  rectangular  matrix.  Multiplying  equation  (5-16)  on  the  left  by 

B  ,  we  nave  Br±S-BTBM,  +  B‘ v 

(5-22) 

After  transposition,  we  get 

BrfiAA, »  B1  (  AS  — 1 ■  v) 

Multiplying  (5-23)  on  the  left  by  (BTB)~'1',  we  have 


AA.-CS'flr'B’CAS-v) 


(5-24) 


3^ 


» 


This  is  the  basic  equation  of  the  best  estimated  value  when  we 

r 

calculate  AX~  with  a  large  amount  of  data.  Statisticians  call  (5-2*0 
the  normalized  or  normal  equation,  and  also  the  wave  filter  equation. 
It  will  reduce  the  effect  of  the  random  error  v  on  If  the  square 

difference  of  each  observation  point  i  of  the  measurement  system  is 
known,  better  results  will  be  obtained  by  weighting  each  point  in 
(5-2*4)  according  to  the  square  difference. 

After  expansion,  it  may  be  shown  that  (5-2*4)  may  be  expressed 
as  the  sum  of  the  corresponding  form  s  of  the  sub-matrices  at  a  point 

^  .  V  x  .V 

AA,-(Brflr'(.£:[AS-v2)=  (£BTS,)“  (£  Af[AS— f‘]  ) 

(5-25) 

Results  similar  to  equation  (5-24)  and  (5-25)  may  be  obtained  by 
using  the  least  square  method  in  detailed  derivations. 

VI.  TOTAL  SYSTEM  ERROR 


If  the  time  that  the  command  for  shutting  off  tne  engine  by  the 

♦» 

guidance  system,  t,  is  before  the  time  t "  when  the  vehicle  thrust 

ends,  then  it  is  necessary  to  delete  th-^  *ffe?t  AX,  of  the  thrust 

f  * 

in  the  interval  (t^,t  )  (i.e.,  the  usually  mentioned  after-effect 
thrust)  on  the  orbital  parameters  from  A>‘  by  ..sing  th--  remotely 
measured  or  externally  measured  parameters  so  "r.a  we  -nay  obtain  the 
orbital  parameter  deviation  caused  by  the  re~al  system  error: 

AA*  =  A>y  —  A/.»  (6-1) 

AX1  is  replaced  by  the  AX,,  calculated  with  '5-25'.  The  total  guid¬ 
ance  error  or  synthetic  guidance  error  is 


(5-1) 

?he  total  guid- 


•  -  1 


(b-2  » 


AQ  denotes  the  circular  error  or  the  deviation  of  the  transverse  and 
longitudinal  impact  points.  The  partial  derivative  is  found  from 
the  calculation  of  deviation  orbit. 
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VII.  CONCLUSIONS 

The  main  ideas  of  the  "ba^ic  principles",  equations  (4-6), 

(6-1  and  (6-2)  and  the  measurement  requirements  in  this  paper  had 
been  published  in  1975.  In  recent  years,  large  volumes  of  experi¬ 
ments,  the  measurement  of  orbital  parameters  during  propulsion-f ree 
stage  and  calculation  of  the  impact  point  deviation  indicate  that 
the  above  theory  is  correct,  and  may  be  used  as  an  auxiliary  method 
to  judge  the  accuracy  of  the  guidance  system. 

From  the  diagram,  it  can  be  obviously  seen  that  the  data  for 
the  descending  stage  will  yield  better  analytical  results  than  those 
for  the  ascending  stage. 

The  proposal,  organization,  implementation  and  perfection  of 
this  method  is  the  joint  effort  of  many  comrades,  including  Wang 
Meizhi,  Wu  iiren,  Zhou  Taiping,  etc.  The  author  only  made  some 
derivations  and  explanations. 
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COMPUTER  SIMULATION  OF  LIQUID 
ROCKET  ENGINE  TRANSIENTS 


Wang  Kechang  03^ 

Abstract 

In  this  papier  the  transients  are  simulated  by  means  of  the  characteristic  method  and 
the  numerical  computation  method  for  start  and  shutdown  processes,  pulse  operation,  four- 
engine  parallel  operation  and  failure  state  of  liquid  rocket  engine  Simulation  is  made  on 
the  computer  ‘‘441B-lir\  The  simulated  results  show  that  in  the  engine  development  com¬ 
puter  simulation  would  be  necessary. 


TA3LE  OF  SYM30LS 

A  channel  cross-sectional  area 

a  sound  velocity  of  propellant 
in  the  channel 

c^,c2,c^  coefficients 

D  channel  diameter 

E  elastic  modulus 

f  frictional  coefficient  of 
flow 

g  gravitational  acceleration 
K  fluid  bulk  modulus 


0  volume  flow  rate 
?  pressure 
t  time 

e  channel  wall  thickness  density 
r  specific  gravity 
3  channel  inclination  angle 

p  density 


I .  INTRODUCTION 


the  operational  status  of  liquid  propellant  rocke~  no  1  r. e c  c" 
various  types  consists  of  pulsed  operation,  multiple  start,  multi¬ 
ple  engine  parallel  operation,  variable  propulsion  operation,  etc. 
For  various  types  of  engines,  it  is  required  that  they  have  good 
dynam_c  properties  during  start,  shut  down  and  propulsion  adjustmen 
which  is  to  say  that  they  should  have  satisfactory  transient  opera¬ 
tional  processes.  Hence,  the  transient  processes  of  a  liquid  pro¬ 
pellant  rocket  engine  represent  a  very  important  aspect  of  rocket 
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engine  research.  Many  people  have  long  been  searching  for  the  laws 
of  nature  in  this  area.  However,  because  the  study  of  these  trans¬ 
ient  processes  involves  the  dynamic  equations  of  the  various  parts 
cf  the  engine,  where  the  channel  equation  of  the  supply  system  is  a 
system  of  hyperbolic  partial  differential  equations,  it  was  very 
difficult  to  carry  out  such  kinds  of  computational  research  work 
before  the  wide  spread  use  of  the  electronic  computers.  It  was  in 
the  60's,  following  the  appearance  of  the  electronic  computers,  that 
ground  was  broken  in  implementing  the  kind  of  calculations  mentioned 
above.  Many  papers  describing  computational  study  in  this  area 
appeared  between  1968-1970  with  satisfactory  results. 


R.  J.  Dorsch  and  D.  J.  Wood  [1]  used  the  wave-plane  method  to 
analyze  the  nonstaticnary  channel  flow  and  computed  the  parameters 
of  the  propellant  supply  system  of  the  liquid  propellant  rocket 
engine . 


John  J.  Boehnlein  [2] 
the  channel  wave  equation, 
channel  and  computed  the  st 
pellar.t  rocket  engine  with 


et  al  used  the  Laplace  transform  to  solve 
analyzed  the  nor.statior.ary  flow  in  the 

u  r* i. t  c  r 00s  3  o*''  t iYip  ]_*  c  v>n— 

satisfactory  result. 


J.  C.  Eschweiler,  H.  W.  Wallace  [3],  ?•  Thompson,  T.  J.  Walsh 
[9]  used  the  characteristics  linear  method  to  solve  the  channel  non¬ 
stationary  flow,  and  simulated  the  transient  operational  process  of 
the  liquid  propellant  rocket  engine  on  the  computer  with  good  results. 


In  tnese  papers,  simple  mathematical  models  are  used  to  treat 
the  combustion  chamber  dynamic  equation.  Hut  W.  T.  Webber,  W.  A. 
Gaubatz  [5]  did  a  lot  of  work  on  the  model  cf  the  combustion  chamber 
while  comrades  Yang  3enlian  and  Chen  Guo^ai  did  a  large  amount  of 
valuable  computational  work  [6]  with  reference  [5]  as  the  foundation. 
In  this  paper,  we  also  start  with  the  characteristics  line  method 
with  a  list  of  the  various  commonly  used  boundary  conditions  in 
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liquid  propellant  rocket  engines , and  simulated  on  the  44lb-I  com¬ 
puter  the  transient  processes  of  start,  shut  down,  pulsed  operation, 
multiple  engine  parallel  operation  and  engine  fault  state,  result¬ 
ing  in  a  list  of  simulation  curves  for  the  various  situations. 

II.  BASIC  EQUATIONS 


The  basic  equations  used  in  analyzing  the  engine  transient  pro¬ 
cesses  are  mainly  the  differential  equations  for  the  nonstationary 
fluid  flow.  [7]  has  given  a  detailed  derivation.  Here  we  shall 
just  list 

0  aQ  .  aQ  A  aP  .  / 


a  a x  +  at 


+  F  fx  +  -AgaaS^O 


\  e 


-1) 


the  continuity  equation: 


L/9.  \ ,  «*  dQ 

p\a  ax  xr-° 


(2-2) 


This  is  a  set  of  hyperbolic  partial  differential  equations  for 
which  it  is  generally  rather  difficult  to  find  accurate  closed  solu¬ 
tions.  [7]  gave  a  detailed  description  of  solving  these  equations 
with  the  characteristics  method.  We  shall  only  list  the  final 
results : 


Pr-P<  +  ~~(Q'-QA)‘  fPa 


+  201^*  ‘  ^  =  0 

A  A  —  a  At 


k:* 


(2-3) 


AA  - ■  aAf  1 


(2-14) 


where  P^,Pg,Pn  represent  respectively  the  pressure  at  points  A,  B 
and  ? . 


QA,Qg,  represent  respectively  the  bulk  flow  rate  of  the 
points  A,  B  and  P  in  Figure  1.  ^pa,P8’^A’®B  are  known •  They  may  be 

determined  from  the  initial  conditions  or  calculated  from  the  pre¬ 
vious  time  interval). 
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It  should  be  noticed  that  the  first  equation  in  Equation  (2-3) 

is  only  true  along  the  +ve  characteristic  line  AP .  Similarly, 

the  first  equation  of  equation  (2-4)  is  only  true  along  the  -ve 

characteristics  BP.  As  there  are  only  two  unknowns  Q  , P  in 

P  P 

equation  (2-3)  and  equation  (2-4),  therefore  and  Pp  may  be  found 
from  equation  (2-3)  and  equation  (2-4).  In  order  to  facilitate 
the  writing  of  the  computer  programs,  we  give  below  a  computational 
equation  for  expressing  the  parameters  of  the  ith  section  channel 
flow.  Substituting  the  symbols,  as  shown  in  Figure  2  into  the  first 
equations  of  equations  (2-3)  and  (2-4),  we  get 


0„-o.5  0,.,-0(,  +  i(Pw-P„  ) 


2  j)  i  !  Q._i !  ■+•  0,. ,  •  Q,.  ■  I ) 


(2-5) 


where 


ap 


P„-o. 5  P,„,+P..tT 


(2-6) 


QPxlP,i - are  the  flow  parameters  of  the  ith  channel  section  at 

the  time  instant  to  be  computed 

0._, .  P,_, - are  the  flow  parameters  of  the  i-lth  channel  section  at 

the  time  instant  previous  to  that  being  computed 
(2-5), (2-6)  are  tj-,e  folow  parameters  of  the  I  +  lth  channel  section 
at  the  time  instant  previous  to  that  being  computed. 
Equations  (2-5),  (2-6)  are  the  computational  equations  that  we  will 
use  frequently  from  now  on  to  compute  the  mid-point  flow  parameters. 
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Figure  1.  Characteristic  grid 


Figure  2.  Grid  for  finding  the 
mid-point 


BOUNDARY  CONDITIONS 


Equations  (2-5),  (2-6)  can  only  be  used  to  calculate  the  flow 
parameters  at  the  mid-point  of  each  computed  channel.  But  we  shall 
discover  that  at  either  end  of  the  computed  channel  there  is  only 
one  characteristics  equation  (C  or  C~ )  that  can  be  used.  The  C~ 
characteristics  equation  can  only  be  used  upstream  while  the  C 
characteristics  equation  can  only  be  used  downstream.  However,  with 
one  equation  it  is  impossible  to  find  two  unknowns  Thus, 

it  is  necessary  to  establish  an  equation  to  connect  the  two  unknowns. 
The  boundary  condition  problem  is,  therefore,  involved.  The  setup 
of  the  boundary  conditions  and  flow  system  Is  closely  related  to  the 
forms  of  the  object.  Detailed  descriptions  of  various  boundary 
conditions  for  the  general  hydroelectric  power  systems  have  been 
presented  in  [3]  ar.d  [91-  But  as  shown  in  Figure  3,  the  liquid  pro¬ 
pellant  rocket  engine  system  is  a  complicated  system  with  boundary 
conditions  of  diverse  types. 

References  [7,11,12]  have  derived  in  detail  the  various  bound¬ 
ary  conditions  of  liquid  propellant  rocket  engine  s-otere  .  These 
boundary  conditions  involve  the  following  principal  objects.  They 
are:  the  propellant  storage  tank,  the  valve  or  flow  regulating 

orifices  in  the  channel,  the  sprayer,  the  propellant  liquid  collect¬ 
ion  chamber,  channel  connectors  (3  channel,  5-channel,  etc.),  cooling 
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jacket,  energy  storage  unit,  centrifuge,  etc.  We  shall  only  point 
out  a  few  important  points  here: 

1.  All  boundary  conditions  are  treated  similarly.  They  are  mainly 

applied  to  the  C+  or  C~  characteristics  equations.  Then  the  equation 

connecting  Q  ,P  ,  (such  as  the  sprayer  equation,  the  valve  equation 
P  P 

or  the  centrifuge  equation,  etc.)  is  listed  based  on  the  shape  of 
the  object  at  the  boundary.  The  equations  are  then  solved  simultan¬ 
eously  by  using  the  continuity  relations  to  obtain  the  pressure  (Pn) 

and  flow  rate  (Q  )  at  the  boundary. 

P 

2.  The  equation  thus  solved  is 
usually  a  quadratic  equation  with 
one  unknown.  The  solution  can  be 
easily  found.  Only  the  boundary 
condition  of  the  energy  storage 
unit  gives  a  nonlinear  equation 
which  we  solved  with  the  Newton 
approximation  method. 


3.  The  boundary  condition  equa¬ 
tions  listed  in  [7,11,12]  have 
been  used  in  extensive  computa¬ 
tions  on  the  computer.  The 
results  show  that  these  equations 
can  be  used. 

IV.  COMPUTER  SIMULATION  OF  VARIOUS  TRANSIENTS  OF  A 
LIQUID  PROPELLANT  ROCKET  ENGINE 

1.  Determination  of  time  increment  At  and  number  of  section  N. 

The  supply  system  of  the  liquid  propellant  rocket  engine  is  a 
complicated  system,  composed  of  many  channels  and  the  lengths  of  each 
channel  may  differ  greatly.  Before  proceeding  with  the  simulation, 
the  problem  of  determining  and  selecting  the  time  increment  At  and 
how  to  divide  the  channels  into  sections  is  encountered. 


Figure  3-  Simplified  diagram 
of  the  supply  system  of  two 
elements  liquid  fuel  rocket 
engine . 


Firstly,  the  time  increment  for  each  channel  should  be  the 
same  so  that  the  boundary  conditions  can  be  used.  At  the  same  time, 
the  Courant  stability  condition  [9]  must  be  satisfied  in  selecting 
the  time  increment,  namely: 

_A/_  )  (^-1) 

Ajc  ^  a 


This  indicates  that  Ax  must  be  greater  than  or  equal  to  a. At. 
Hence,  the  selection  of  At  under  the  condition  that  the  speed  of 
sound  is  already  determined  is  mainly  limited  by  Ax  while  Ax  must 
be  constrained  by  the  shortest  channel  in  the  system.  If  the  short¬ 
est  channel  is  short,  then  Ax  will  be  small  and  At  cannot  be  taken 
very  large.  When  At  is  small,  the  computational  time  and  expense 
will  increase  so  that  for  a  large  number  of  channels,  the  fine  divi¬ 
sion  nay  even  exceed  the  memory  of  the  computer.  Based  on  our  exper¬ 
ience,  we  think  that  At  should  be  chosen  as  follows.  Firstly,  the 
system  channels  should  be  analyzed  to  find  the  shortest  one.  If 
this  channel  is  indeed  too  short,  then  it  should  be  treated  separately 
by  the  finite  difference  method  similarly  to  the  procedure  for  the 
sprayer  liquid  collector  chamber.  If  this  short  channel  is  not  too 
short,  then  it  may  be  treated  in  two  sections,  from  which  a  At  value 
will  be  found. 

=  (9-2) 

This  At  then  can  be  used  for  the  whole  system. 

However,  we  must  notice  that: 


1)  After  At  is  determined  by  analyzing  the  shortest  channel,  in 
general,  it  will  not  change  any  more.  Then  the  other  channels  may 
be  treated  by  the  following  method: 
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N  is  the  number  of  divisions  of  the  ith  channel  and  should  be 
an  integer.  But  the  N  value  found  from  equation  (9-3)  is  often  non¬ 
integer.  Hence,  both  V.  L.  Streerer  and  M.  H.  Chaudry  have  pointed 
out  that  due  to  the  approximate  nature  in  the  calculation  of  a,  i.e.. 
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its  value  cannot  be  accurately  known,  therefore,  it  is  possible  to 
adjust  the  value  of  speed  of  sound  approximately  to  make  N  value  an 
integer . 


2)  We  have  carried  out  comparison  c  • mput at  ions  on  the  computer  on 
the  effect  of  the  number  of  divisions  of  the  shortest  channel  and 
discovered  that  when  N  was  taken  to  be  2  and  10  the  results  were 
very  close,  which  indicates  that  very  accurate  results  are  attain¬ 
able  when  the  number  of  divisions  is  reduced.  Table  1  lists  the  data 
of  the  analytical  calculations  of  engine  disturbance  for  comparison. 


TABLE  1.  Comparison  of  the  calculated  values  of  the  r  rir.  c  ip  a  1 
parameters  5)  for  difference  H 
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io  ms 
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0.9819 
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i 

0  999 
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|  1  0018  | 
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0  9995 
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1  001 

During 
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calculatio 

n,  we  found 

that 

the 

effect 
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the  number 

of  di 

visions  is 

very  attractive 

as 

long 

as 

satisfactory 

accuracy  can  be 

maintained 

When  N  = 

the 

t  ime 

in 

crement  of  the 

disturbance  computations  is  At  =  0._  ms.  To  compute  one  At,  it  would 
take  55  seconds  of  computational  time  on  the  441B-I  computer,  while 
for  N  =  2,  the  time  increment  At  =  0.5  second  and  it  only  took  about 
10  seconds  to  compute  one  At.  The  total  computational  time  for  one 
engine  disturbance  was  also  reduced  from  three  hours  to  half  an  hour. 
Thus,  the  choice  of  N  (hence  the  determination  of  At)  is  a  very 
important  question,  worthy  of  our  attention. 

2.  Simulation  of  liquid  propellant  rocket  engine  start  and  shut-down 

Start  and  shut-down  are  two  important  transient  processes  of  a 
liquid  propellant  rocket  engine.  They  are  related  to  the  reliability 
of  start,  reliability  of  stage  separation,  etc.  We  have  carried  out 
calculations  on  the  441B-I  computer  :•  f  various  types  of  start  an':  oh  ; 


Figure  u .  curve  of  combustion  chamber  pressure  (P  )  with 
time  with  the  valve  right  next  to  the  injector  andCno 
.Liquid  collecting  chamber 
-no  ignition  delay;  2 — no  ignition  delay 


valves  versus'time1"  ^  ?reSSUre  in  fr°nt  °f  the  e™?*11"* 
pressure  in  front  of  the  oxidizer  valve; 
pressure  in  front  of  the  combustible  valve 
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down  processes  for  an  attitude  control  engine  v-ith  11.35  kg  of 
thrust,  using  +  R-NK-NH-50  as  propellant.  The  results  of  the 

computation  are  shown  as  the  curves  in  the  following  diagrams. 

Analysis  on  the  result  may  be  found  in  [?]. 

3.  Simulation  of  the  attitude  control  engine  pulsed  operation  and 
4  engine  parallel  operation  transients. 

The  attitude  control  engine  used  on  flying  vehicles  has  the 
following  characteristics.  Firstly,  the  number  of  attitude  control 
engines  is  large,  usually  4  or  more,  and  they  all  use  the  same  pro¬ 
pellant  supply  channel.  As  the  attitude  control  engines  do  not 
usually  start  or  shut-down  at  the  same  time,  therefore,  one  must 
clearly  understand  the  effect  of  start  or  shut-down  of  one  or  two 
engines  on  the  operation  of  the  other  engines.  The  other  character¬ 
istics  of  attitude  control  engines  is  the  large  number  of  re-starts. 

For  instance,  the  number  of  starts  of  the  attitude  control  engines 
on  the  Apollo  may  reach  100,000  and  the  operation  frequency  4q  times/ 
sec.  What  is  the  effect  of  adjacent  pulse  operations  on  one  another 
under  this  kind  of  operation  frequency?  In  order  to  understand  this 
problem,  we  separately  simulated  the  following  cases: 

1)  Restart  10  ms  after  the  end  of  the  first  pulse.  Here  the  pressure 
in  the  oxidizer  is  in  the  low  pressure  wave  region  and  that  in  the 
combustible  channel  is  in  the  high  pressure  wave  region. 

2)  Restart  30  ms  after  the  end  of  the  first  pulse.  Here  the  pressures 
in  both  the  oxidizer  and  combustible  channels  are  in  the  low  pressure 
wave  region. 

3)  Restart  60  ms  after  the  end  of  the  first  pulse.  Here  the  pressure 
in  the  oxidizer  channel  is  Ln  the  high  pressure  wave  region  and  that 
in  the  combustible  channel  is  in  the  low  pressure  wave  region. 

4)  Restart  100  ras  after  the  end  of  the  first  pulse.  Here  the 
pressures  in  both  the  oxidizer  and  combustible  channels  are  in  the 
high  pressure  wave  region. 
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5)  Effect  of  channels  of  different  diameters  on  the  start  process. 

6)  Two  engines  start  first  and  after  several  ms,  the  other  two 
engines  will  then  stare. 

The  results  of  the  simulation  are  shown  in  the  following  curves. 
Analysis  of  the  results  may  be  found  in  [11]. 

A.  Simulation  of  liquid  propellant  rocket  engine  fault  state  transient 

Analysis  of  liquid  propellant  rocket  engine  fault  state  is  very 
important  and  very  complicated.  Fault  analysis  is  involved  from  the 
beginning  to  the  end  in  the  test  analysis  of  the  engine  on  test 
stand  and  in  the  flight  test  analysis  of  the  engine.  In  the  test 
stand  testing  stage,  it  is  still  possible  t:  use  th-  high  speed  photo¬ 
graphy,  experimental  parameter  curves  recorded  during  the  test  and 
the  engine  itself  after  the  test  for  a  synthetic  analysis.  In  the 
flight  test,  because  it  is  difficult  to  find  complete  remains  of  the 
engine,  we  can  only  make  cur  analysis  based  on  the  remote  measurement 
data.  Hence,  if  one  can  simulate  with  computation  the  transient 
processes  of  various  important  parameters  of  the  engine  (e.g.,  com- 

^  ~  Z  n  0  2LT  *0  0  1"*  C  v‘» <=*  f  3  't  ^  r*  GT^  0  ^  2.  n  *"  °  **  V1  c.  ^  0  'G  v*>  L,  1.  P  ^  v-1 p  m  - 

speed,  etc.)  under  various  fault  conditions  (e.g.,  channel  leak, 
injector  blockage,  etc.),  the  analysis  and  research  on  engine  fault 
will  be  greatly  benefited  since  we  only  need  to  compare  the  transient 
curves  of  various  parameters  actually  recorded  during  engine  fault  to 

r'  T  c-  ',{'*■  r.  1  ^  -5  ^  u*"  '  vT  "  '  a  *■)'  f  ^  ^  ^  •  •  ~4  v-»  <  ,  o  -+  p  "  '"0  j  ^  I  *" 

conditions  in  order  to  determine  the  type  and  location  of  the  fault. 
This  will  help  in  narrowing  down  the  scope  of  analysis  and  reducing 
the  work  load  for  analysis,  and  hence  will  promote  greatly  the  pro¬ 
gress  of  engine  development. 


The  computation  of  engine  fault  transients  is  primarily  dependent 
upon  the  establishment  of  fault  state  computational  model.  In  [12], 
an  engine  fault  computational  model  of  a  medium  thrust,  pumped 
propellant  with  2  engine  group  elements  is  presented.  The  following 
diagrams  show  the  computational  results.  An  analysis  of  the  results 
may  be  found  in  [12]. 
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1 — pump 
power; 

6  — fuel 


gure  12.  Transient 
rotational  speed; 

4 — fuel  pump  power; 


during  oxidizer  channel  leakage 
2 — turbine  efficiency;  3 — oxidizer  pump 
5 — combustion  chamber  fuel  flow  rate; 


flow  rate  in  fuel  vaporizer 
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Figure  13  •  Parameter  variation  during  oxidizer  valve 
1 — fuel  vaporizer  pressure;  2 — combustion  chamber  pressure; 

3 — oxidizer  flow  rate  in  fuel  vaporizer;  4 — combustion  chamber 
oxidizer 

V.  CONCLUSIONS  AND  PROPOSALS 

1.  Computer  modeling  of  liquid  propellant  rocket  engines  is  very 
instructive  in  the  development  of  the  engine.  It  helps  us  in  select¬ 
ing  the  principal  parameters,  in  proving  the  major  plans  and  in  pre¬ 
dicting  the  whole  transient  process.  Today  when  the  application  of 
the  electronic  computer  is  ever  widening,  computer  modeling  will 
become  an  Important  link  in  the  design,  development  and  experiment¬ 
ation  of  liquid  propellant  rocket  engines. 

2.  The  keys  to  the  computer  modeling  of  engine  transients  are  (1) 
the  establishment  and  gradual  perfection  of  the  mathematical  model 
of  various  parts  of  the  engine.  For  instance,  there  are  now  more 
than  a  dozen  computational  models  for  the  combustion  chamber  computa¬ 
tional  model  internationally.  We  should  study  and  analyze  these 
models  to  establish  a  mathematical  model  that  fits  our  needs.  (2) 
Development  of  the  commonly  used  computer  programs  (including  the 
study  of  the  computational  methods).  (3)  Supply  of  experimental 
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data  for  the  various  parts.  For  instance,  it  is  still  difficult  to 
find  the  flow  coefficient-time  curve  and  the  degree  of  openness  vs. 
time  curve  during  valve  opening  and  closing.  Much  work  needs  to  be 
done  before  computer  modeling  of  engine  transients  can  be  perfected 
and  used  in  practice.  In  this  paper,  we  have  only  provided  a  rough 
outline.  We  welcome  your  criticism  so  that  the  worK  of  computer 
modeling  may  be  further  developed. 


REFERENCES 

( 1  ]  R.G.  Dorset,  DJ.  Wood,  Distributed  Parameter  Analysis  of  Pressure  and  Fo*  Disturbances  in  Rocket 
propellant  feed  systems,  NASA  TN  D-3529,  1966. 

(2)  John  J.  Boebnlein,  ct  al,  Generlized  Propulsion  System  Model  for  NASA  Manned  Spacecraft  Center, 
NASA-CR-1 14915,  1971. 

[  3  ]  J.C.  Eschweiler,  H.W,  Wallace,  “Liquid  Rocket  Engine  feed  System  Dynamics  by  the  Method  ofCharac- 
terisbea,'’  Transactions  of  the  ASNE  Series  B  Vol.  90.  NO  4.  1968 

(4]  P.F,  Thompson,  T.J.  Welsh,  Characterization  of  Attitude  Control  Propulsion  Systems.  NASA-CR- 
115183.  1971. 

(5]  W.T.  Webber, kW. A.  Gaubatz,  “Calculation  of  the  Ignition  and  Start  Transient  of  Liquid  PropeOant 
Rocket  Engines,"  WSS/d  paper  70-23. 

[6]  Xiang  Ben-lian,  Chen  Guo-tai ,  Computer  Calculation  cf  Hot 
Start  and  Shut  Down  of  Liquid  Fuel  Rocket  Engine,  Journal 
of  Engineering,  1979,  a 

[7]  Wang  Ke-ehang,  Computation  of  the  Start-up  and  Shut-down 

of  Liauid  Fuel  Rocket  Engine,  Journal  of  Engineering,  1979,  9 

L8J  E.B  Wylie.  V.L.  Streeter,  Fluid  Transients.  1978. 

(24  M  H.  Chauanry,  Applied  Hydraulic  Trurisionts.  I9'9. 

(101  D.J.  Wood,  S.E.  Jone.  “Water  Tammer  Charts  for  Various  Type  of  valve*  ASCE  Vol  99  HY  I  m 
167178.1973.  ’  ,p|> 

[11]  Wang  Ke-chang,  Some  Applications  of  Numerical  Methods  in 
Calculating  the  Transients  of  Liquid  Fuel  Rocket  Engines, 
Journal  of  Engineering,  1979,  a 

[12]  Wang  Ke-chang,  Computer  Model  Analysis  of  the  Fault  State 
of  Pump-Pressured  Liquid  Fuel  Rocket  Engine 


53 


» 


REPEATED  TRIAL  OF  GUIDANCE  IN  SPACE  VEHICLES 
He  Lucheng 


ABSTRACT 

This  article  covers  the  advantages  of  the  repeated 
method  and  the  applications  of  the  repeated  method  in 
strap-down  inertial  system.  The  differences  of  this 
method  with  those  described  in  the  previous  article  [1], 
[2]  are: 

1.  Not  only  the  reliability  of  the  system,  but  also 
the  accuracy  of  the  system  will  be  improved. 

2.  The  key  components  which  affect  the  accuracy 
in  guidance  should  be  repeatedly  fit  in  the  optimum 
directions  and  those  which  affect  reliability  should  be 
considered  to  be  in  the  orthogonal  or  non-orthogonal 
repeated  method. 

3.  This  article  provides  the  method  of  inspecting 
failure  by  means  of  range  and  derives  the  formulae  of 
calculating  accuracy  and  effectiveness.  This  method 
will  make  calculation  simple  and  work  of  the  system  reli- 
ab  le . 


^ .  Applying  the  theory  mentioned  above,  the  error 
in  guidance  will  be  decreased  by  and  the  ineffect¬ 

iveness  can  also  be  decreased  by  two  orders  of  magnitude, 
if  a  longitudinal  accelerator  and  a  two-freedom  gyro  are 
added  to  the  strap-down  guidance  device. 

I .  Formulation  of  the  problem 

Repeated  trial  methods  may  be  used  to  enhance  the  accuracy  of 
space  vehicle  guidance  by  compensating  for  systematic  errors  and 


reducing  random  errors.  The  repeated  trial  method  is  not  simply  to 
enhance  the  accuracy  by  the  a,  method.  It  also  comDensates  for 
the  unknown  systematic  errors  (gyro  acceleration  meters  are  in¬ 
stalled  with  one  normal  and  the  other  reversed  and  backward,  so 
that  when  the  data  from  the  two  meters  are  averaged,  systematic 
errors  are  compensated)  and  cancels  rolling  error.  The  more  import¬ 
ant  advantage  of  the  repeated  trial  method  is  to  improve  reliability, 
as  unreliability  is  potentially  more  harmful.  The  strap-down  method 
eliminates  the  constant  level  platform  and  its  stabilizing  system. 

As  electronic  computers  are  being  perfected  and  miniaturized,  the 
weight  and  size  of  the  strap-down  inertial  system  may  be  made  very 
small.  The  most  important  advantage  of  the  strap-down  inertial 
system  is  the  ease  in  using  the  repeated  trial  method  which  greatly 
improves  reliability. 

II.  Several  reoeated  trials  method 


1.  Arithmetic  average  method 

In  order  to  understand  the  advantages  and  disadvantages  of 
various  repeated  trial  methods,  we  shall  first  describe  the  error 
probability  model  of  a  single  inertial  assembly.  For  example,  the 
error  model  of  a  gyro  is 

5=  V («.,  +  «„)!.",  =  3,.  if  4-*, W  (  I  ) 

•m'  * 

where  means  vector  dot  product,  i.e.,  ^  3 

_  •  ■  1 

w,  is  the  environmental  factor,  including  the  overload  acceleration 
experienced  by  the  sensor,  as  well  as  the  coupling  cf  accelerations 
along  other  directions,  and  such  environmental  stresses  as  tempera¬ 
ture  rise,  vibrations  and  gas  pressure. 

are  the  uncorrected  error  coefficient  values  correspond¬ 
ing  to  the  environmental  stress.  *•  is  the  mathematical  expect¬ 
ation,  referred  to  below  as  regular  error;  5,  is  the  random  part, 
is  the  system  output  error. 
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Under 

ary  to  take 
formed  into 
independent 
error  model 


the  action  of  some  environmental  stresses,  it  is  necess- 
the  absolute  value  of  some  errors  which  can  be  trans¬ 
regular  errors  and  random  errors.  If  the  errors  are 
and  are  composed  of  many  small  disturbances,  then  the 
may  be  written  as 


m=»£[S]  =  Wf-  W  C  •?  ) 

«•! 


Ac 

N(m,  a) 
element 
ground . 


cording  to  the  Central  Limit  Theorem,  normal  distribution 
must  be  satisfied.  Many  trials  are  performed  on  the  same 
under  the  influence  of  each  environmental  stress  on  the 
Each  environmental  sti'ess  is 

fn,  =  r,cr, 


m,  sometimes  may  be  several  times  that  of  ck  . 


We  try  to  find  m  to  make  corrections  by  using  as  much  as  poss¬ 
ible  the  data  of  many  trials  and  from  ground  signals,  but  sometimes 
it  is  not  easy.  However,  if  we  use  an  even  number  of  basically 
identical  elements,  with  half  of  them  normal  and  the  other  half 
reversed  and  backward,  and  then  average  the  measured  data,  then  the 
regular  errors  will  be  mostly  cancelled  and  random  error  reduced.  If 
the  number  of  repetitions  is  n,  then  the  random  error  will  be 
reduced  to  --i,  .  When  this  simple  arithmetic  averaging  method  is 

v  n 

used  on  a  certain  object  without  too  many  measurements,  any  anoma¬ 
lous  data  can  generally  be  detected  by  inspection.  But  for  unmanned 
automatic  guidance  systems,  any  faulty  data  or  data  with  large 
deviations  must  be  instantly  detected.  The  normal  data  are  then 
averaged  so  that  the  accuracy  and  reliability  of  the  system  may  be 
guaranteed.  Otherwise,  the  average  will  be  in  error  even  if  there 
should  only  be  one  piece  of  data  with  fault  or  with  large  deviations 
and  the  whole  system  may  become  ineffective.  In  this  way,  the 
ineffectiveness  may  be  increased  n  times.  The  improvement  of  accuracy 
actually  means  the  improvement  of  the  probability  of  a  successful 
mission,  but  arithmetic  averaging  may  increase  ineffectiveness. 
Sometimes  the  loss  may  be  more  than  the  gain. 
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In  order  to  overcome  the  shortcomings  of  the  arithmetic  aver¬ 
aging  method,  the  0-1  checking  method  is  proposed,  in  which  data 
with  large  deviations  are  automatically  detected  and  isolated. 

The  basic  logic  is  to  permutate  the  data  according  to  software 
programs  and  to  subtract  in  pairs.  When  the  difference  exceeds  a 
certain  value,  1  is  assigned;  otherwise  0  is  assigned.  Then  the 
location  of  an  error  is  found  based  on  the  way  l’s  appear  and  is 
iisciPd©-'!  • 


In  the  strap-down  guidance  system,  an  orthogonal  coordinate 
system  is  fixed  on  the  vehicle.  Each  axis  must  complete  automatic 
fault  checking  with  respect  to  three  fixed  acceleration  meters  and 
three  single-freedom  gyros.  In  this  way,  the  number  of  elements  is 
increased  three-fold.  To  reduce  the  number  of  elements,  a  regular 
dodecahedron  scheme  was  proposed  in  reference  [1]  and  a  non-ortno- 
gonal  system  was  proposed  in  reference  [2].  Both  of  these  will 
reduce  the  number  of  elements  from  three  times  to  two  times  or  less. 

2.  Dodecahedron  method 

Let  Oxyz  be  a  coordinate  system  fixed  to  the  vehicle.  The 

meter  input  axis  is  fixed  to  six  axes  . u,  .  u,  and  u,  are  in 

the  xy  plane  making  angles  -fa  and  -a  with  the  x-axis. 
are  similar  as  shown  in  the  diagram. 


Let  c—cosa,*- sina, 

of  the  vehicle  may  be 


.  Then  the  visual 
analyzed  as 

™  caM—  sar  -f.  l 

“;  =  ca,  say 

I 

B,™ca, —saM  +  bj 

,u,-ca,  +  sa,  +  b, 

j«,-cof-*a.  +  6, 
i»i«ca,  +  sa,  +  6. 


acceleration 


(S) 
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or  in  matrix  notation  h=Ax-t~o 

where  b  is  the  error  matrix,  u  is  the  measurement  matrix,  A  is  the 
transformation  matrix,  x  is  the  ballistic  motion  state.  If  the  u 
data  of  the  various  meters  are 
all  normal,  then  by  the  least 
square  method 

i~(ArA)~'ATi  (&) 

Theoretically,  to  solve  three 

unknowns  with  six  equations, 

a  .a  ,a  ,  a  fault  may  be  per- 
x  y  z 

mitted  to  occur  for  three  mea¬ 
sured  data.  Reference  [1] 
listed  the  fault  data  detection 
method.  For  example,  from  the 
unand  u2data  in  equation  (5), 
we  can  estimate 

2  co*  a  ' 

and  from  u^  and  Ug  in  equation 

3.=  we  can  estimate 

_  M.-M. 

2Aina  * 

If  the  difference  between  the  2  exceeds  the  allowed  range,  then 

one  or  mere  of  u, ,  u„,  uc  and  u^  are  faulty.  The  same  holds  for  the 
l  u  o  o 

other  u’s.  Finally,  the  fault  location  may  be  determined.  If  no 

more  than  three  pieces  of  data  are  faulty  in  equation  (5),  then  the 

computer  program,  will  check  out  the  faulty  equaton  and  discard  it. 

Then  the  least  square  method  is  used  again  to  solve  for  a  ,a  ,a  . 

x  y  z 

3.  Non-orti.igonal  system  method  [2] 

An  oblique  set  of  axes  is  selected  with  unit  vectors  e? ,  e~,  e® 

1  ’  2  ’  3 

which  satisfy 

j  when  i  /  j 

l  . 

i.e.,  the  three  axes  are  neither  colinear  nor  mutually  perpendicular 
This  is  called  a  non-orthogonal  system. 


58 


A  unique  reciprocal  system  c!.c;.c;,  nay  be  determined  for  this 
non-orthogonal  system  satisfying  the  following  relations: 


eo*  a., 
0 


(7) 


when  i  =  j 
when  i/j 

i.e.,  c^  is  perpendicular  to  the  plane  formed  by  e|,  and  form  an 
angle  a  with  e°.  e°  is  perpendicular  to  the  plane  formed  by  c°  and 
c°.  The  other  cases  follow  similarly.  cos  a.,  is  the  direction 
cosine  of  the  corresponding  axis . 


The  gyro  or  acceleration  meter  is  fixed  on  the  vehicle  along 
the  six  axes  «f.«?.ef,c? .cj.c*  .  The  data  measured  are 
respectively.  Assume  that  there  exists  a  space  vector  2.  then  the 
projections  of  the  vector  Hi  along  the  six  axes  are  £„=«*•*,-,  Qt,=*S'C? 
measured  by  the  meters  on  them.  2  may  be  decomposed  into  three 
components  along  «f.*f.«;  directions  with  size  *i.*«.*i  .  Similarly, 

it  may  be  decomposed  along  the  cT.eJ.c*,  directions  with  sizes 

c..c,.c,.  ,  From  equation  (7),  we  know  that  c^  and  c^  are  normal  to 


Hence,  Q.l"c,ca»al .  Similarly 


tu 


cot  a,  0 

0  cot  a. 

Q~  J  L  o  o 


o 

0 

cot  a. 


J 


c, 

C. 

C.  J 


(8) 


Equation  (8)  may  be  written  in  matrix  form 

or  C  **  D~'  Q, 

Similarly  Q,  —  DE ,  or 

E=D~'Q, 


Q.  —  DC, 


(9) 

(10) 


The  angles  aT2»a23  a3l  between  the  three  axes  cf.e5.e5  of  the 
non-orthogonal  system  are  known  and  should  be  original  parameters. 
The  measured  valuesQ,,.  H.-.L),,.  may  also  be  expressed  as  the  sum  of 
the  projections  along  this  measurement  axis  of  the  three  components 
e  .,e  2 ,  s  ^  • 

1  cm  a,, 

co»a,t  1 
cm  a,,  CMS,, 


r* 

-| 

.a.. 

L 

CMC., 

CM,, 

1 


(ID 
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I 


l 


I 


Substituting  equation  (10)  into  equation  (11),  we  get  the  value 

from  which  fi  is  estimated  with  the  measured  value  fi  . 

e  c 


*&,  ■ 

■a, " 

=  TD~' 

a.. 

!  &. . 

.  flrl  - 

.!J  =  DE,  according  to  equation  ( 11 )  ,  E  =  T~'Qe.  hence 


■a,  ■ 

a,  ] 

A. 

- DT~ ' 

a,  | 

_a.  J 

After  the  angles  a^.  of  the  ncn-orthogonal  system  are  determined, 
a^jEpja^  may  be  uniquely  determined  through  solid  geometrical  rela¬ 
tions.  The  difference  between  the  estimated  value  and  the 
measured  value  ft.,  may  be  formed  from  equations  (12''  and  (13) 


£  I  “  £ii — ,  i =  e, .  e, ,  e, .  c,  %cI%c, , 

A  critical  region  >/  may  be  set  uo  so  that  if 

a. 

normal  and  *.=  o  but  if  \*A>Wt  then  it  is 
If  ft  ,  is  faulty,  then  according  to  equations 

C  x 


(U) 

i«.| ,  it  is 
faulty  and  fi”1* 

(12)  and  (13) 


re,.c,.e,.e,  =  i 

«»*«,  -»  0 


The  other  cases  are  similar.  Based  on  the  regularity  of  the 
appearance  of  0-1,  we  can  determine  and  isolate  the  fault  of  any 
gyro. 


The  common  weakness  of  the  dodecahedron  method  and  the  non- 
orthogonal  method  is: 

(1)  The  principal  environmental  stresses  have  not  been  taken 
into  consideration.  For  example,  the  visual  acceleration  along 
the  longitudinal  axis  is  more  important  for  the  vehicle,  while  the 
visual  accerlations  W  Wz^  along  the  transverse  axes  are  very 
small . 

Let  the  constant  drift  of  a  certain  gyro  with  two  degrees  of 
freedom  to  be  l°/hr,  and  its  overload  drift  to  be  l°/hr  g  .  When 
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r 


applied  to  an  I CBM,  the  spread  caused  is  shown  in  Table  1. 


i  *•  j  * 

;  *» 

z  }  K.+K.#r,+Kjt'., 

K.+  K,*„  +  K,W.> 

Km^KxW  gi  +  KtfT  || 

"n 

r  ?!£**«£** 

*• 

K. 

ft 

“  ai 

44 

* 

(£«> 

St 

39.3 

n  • 

3.4 

1.4 

13.4 

1 — direction  of  the  angular  momentum  H;  2 — angular  drift  formula; 

3 — measured  angular  accelertion;  U — the  principal  term  causing  the 
drift;  5 — induced  spread 

From  the  above  table,  it  can  be  seen  that  the  gyro  measures 
w„,  and  the  z,  gyro  measures  w  , .  Due  to  the  effect  of  W  , ,  their 
errors  are  twice  that  of  the  gyro.  If  the  direction  of  the  ang¬ 
ular  moment  H  of  the  x,  gyro  deviates  from  the  x,  axis,  then  the 
•  ±  1 

component  W  ^  will  be  involved,  causing  a  larger  error,  while  the 
errors  for  the  y.^  gyro  and  the  z^  gyro  can  only  be  minimized  in  mea¬ 
suring  a  , .  Similar  situations  hold  true  for  the  acceleration  meter 

Xi 

also.  Accordingly,  large  err:rs  will  be  introduced  through  the  use 
of  a  non-orthogonal  coordinate  system  or  dodecahedron  which  should 
only  be  used  when  accuracy  is  not  a  concern.  The  xev  element  affect¬ 
ing  guidance  accuracy  (such  as  longitudinal  acceleration  meter' 
should  have  multiple  installations  along  the  optimal  direction  sc 
that  accuracy  and  reliability  may  be  effectively  enhanced.  The 
transverse  acceleration  meter  is  not  a  key  element  and  may  be  cut 
off  when  it  becomes  faulty  without  much  effect.  Hence,  there  is  nc 
need  for  duplication. 

(2)  Use  of  dodecahedron  or  r.on-orthcc:onal  coordinate  system 
methods  may  easily  produce  a  wrong  diagnosis.  If  two  fault  data 
are  not  very  different  among  n  pieces  of  data  or  one  piece  of  data 
with  large  deviations  is  not  much  different  from  another  piece  of 
normal  data  but  very  different  from  other  pieces  of  data,  then 
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confusion  will  appear  in  the  0-1  combinations.  In  particular,  when 
installation  directions  are  not  the  same,  some  meters  are  more  sen¬ 
sitive  to  environmental  stresses,  thus  causing  even  more  confusion. 

(3)  The  computation  is  more  complicated  when  one  uses  the 
dodecahedron  or  the  non-orthogonal  coordinate  system  method. 

Due  to  the  above  mentioned  causes,  we  propose  in  this  paper 
that  by  installing  multiple  units  of  key  elements  along  optimal 
directions,  both  accuracy  and  reliability  will  be  enhanced.  Atti¬ 
tude  stabilising  systems  may  be  duplicated  orthogonally  or  non- 
orthogonally  to  enhance  only  reliability  without  duplicating  the 
secondary  parts.  The  method  for  fault-checking  will  be  described 
below: 

III .  Fault-checking  by  polar  difference 

1.  Principle  of  checking  by  polar  difference 

There  are  generally  two  forms  of  gyro  error,  one  with  no  output 
(broken  line,  broken  electric  source)  and  the  other  the  fast  rota¬ 
tion  phenomenon  (damaged  meter  or  appearance  of  extremely  large 
positive  or  negative  deviations).  The  first  type  of  faults  is  easy 
to  check.  If  there  is  no  positive  or  negative  pulse  from  a  dynamic 
pressure  gyro  (no  voltage  after  rectification),  then  it  means  a 
broken  line  or  broken  electric  power  source.  To  check  for  the  fast 
rotation  phenomenon,  it  is  simpler  to  use  the  polar  difference 
method . 

Under  the  same  conditions,  because  the  data  have  random  errors, 

the  largest  value  x  „  and  the  smallest  value  x  ,  of  n  measured 
°  max  min 

values  » • • ■ » xn  are  also  random.  The  polar  difference  W  for  this 

set  of  measurements  is 

The  sample  polar  difference  is  a  random  variable.  Let  it  be  after 


after  normalization. 


Let  the  real  values  cf  a  specific  set  of  measurements  be 
represented  by 

=  (*.«-  (15) 

Now  let  a  certain  critical  region  W  be  defined.  If  ^>H». 

3. 

then  we  consider  the  data  to  be  faulty.  If  ,  then  we 

consider  the  n  nieces  of  data  are  all  normal. 


Under  the  condition  that  the  volume  is  n,  the  probability  model 

that  the  random  polar  difference  W  is  greater  than  some  assumed 

o 

value  W  may  be  derived  with  a  multinomial  distribution  formula. 


Let  the  probability  density  function  of  x  be  f(x).  Let  n  data 
be  measured,  any  one  of  which  has  the  same  probability  of  falling 
on  the  smallest  value.  Let  this  value  be  x.  Its  probability  is 
nf(x'.  The  other  n-i  data  will  all  fall  in  the  range  [x,x+W]  with 
probability  Ru)du  }  ’  .  As  x  varies  in  t.ne  range  — oo~  +  oc 

then  the  probability  that  n  data  satisfy  the  polar  difference  is 


__/(*)  I  4  f(u)du 


J 


:t s  complementary  event  is  that  one  or 


more  of  the  n  data  will  fall  outside  the  polar  difference  W  with 


probability  .  Let  it  be  a,  i.e., 

rvrxr. I.)- 1-.  ('_/<»>  [["‘7c«)J.  V 


Oi) 


In  the  above  equation,  if  a  is  specified,  then  there  is  only 

one  unknown  Wa  left.  The  result  of  solving  the  equation  is  shown 

in  Table  2  [  3  j  -  That  is  to  say  that  with  normal  data 

If  H  ,  then  we  have  a  small  probability  event  (one  possibility 

is  sudden  ineffectiveness,  and  the  other  is  large  deviation).  We 

should  then  discard  the  cair  cf  data  x  .  , x  .  One  of  the  two 

min  max 

pieces  cf  data  may  be  a  normal  one.  Wouldn't  it  be  a  waste  to  dis¬ 
card  it?  It  may  be  proved  in  sequential  statistical  theory  that  when 
the  errors  of  the  peripheral  terms  are  large,  there  is  no  advantage 
in  including  them  for  accuracy  calculations.  On  the  other  hand,  by 
discarding  them  we  may  achieve  two  goals:  (1)  reliability  will  be 
improved  by  discarding  ineffective  data;  (2)  accuracy  will  be  enhanced 
by  discarding  data  with  large  errors.  After  the  first  round  of 


TABLE  2 

W,  n 

,  ■  ! 

-»«. . 

*» 

* 

3 

4 

5 

6 

7 

8 

9 

0,05  ' 

} 

2  77 

3  31 

3  63 

3.87 

4.03 

4.17 

4.28 

4  39 

0.10 

2.33 

2  90 

3  34 

3  48 

3  6b 

3.81 

3  92 

4  04 

0  15  I 

2.04 

2.63 

2.98 

3.23 

3  42 

3.57 

3.70 

3  81 

checking,  rhe  ineffectiveness  will  have  been  lowered  by  several 
orders  of  magnitude.  For  if  to  be  further  lowered,  the  remaining 
data  may  be  checked  for  a  second  time  with  critical  region  deter¬ 
mined  by  the  method  in  C  3 1  - 

For  example,  when  sample  volume  n  =  4  (with  four  acceleration 
meters),  take  a  =  0.1  to  get  W  =  3-24  from  the  table.  Thus,  the 

3. 

probability  that  the  polar  difference  is  greater  than  3.24  o  should 
be  less  than  10".  Hence,  we  discard  this  pair  of  data.  If  the 
sensors  possess  a  regular  error  yo  with  half  installed  normally  and 
the  other  half  reversed  and  backward,  then  the  critical  region  may 
be  taken  as  W  =  (3-24  +  2y)a. 

a  is  to  be  determined  according  to  the  actual  situation.  If  a 
is  too  large,  then  type  1  error  will  be  made  (treating  normal  data 
as  abnormal)  with  a  high  probability,  while  if  a  is  too  small,  then 
it  Is  easy  to  introduce  data  with  large  deviations,  therefore,  com¬ 
mitting  a  type  2  error. 

After  polar  difference  checking  with  faulty  data  or  data  with 
large  deviation  discarded,  the  unreliability  defect  of  the  arith¬ 
metic  averaging  method  will  be  solved. 

2.  Ineffectiveness  checking  using  polar  difference 

In  n  pieces  of  data  measured, by  n  instruments,  at  least  (n+l)/2 
(take  integer  value)  will  be  ineffective  and  at  least  (n+l)/2 
ineffective  data  will  all  be  either  greater  or  smaller  than  the  real 
values,  i.e.,  the  system  will  fail  when  all  the  errors  are  one-sided 
(Note:  this  situation  is  hard  to  check  even  with  0-1  checking  since 
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the  data  with  a  large  deviation  may  be  close  to  one  another).  The 
probability  that  each  data  is  larger  or  smaller  is  respectively  1/2 
For  (n+l)/2  or  at  least  (n+l)/2  ineffective  data  to  be  one-sided, 
the  formula  for  the  total  probability  that  the  system  be  ineffect¬ 
iveness  is 


0,= 


(1-9) 


(1-9) 


it 


»-t-l 

2 


]- 


1  -r 


[~“h* 


)  + . +  9’ 


(17) 


Ordinarily  n  is  not  too  large.  When  the  ineffectiveness  of  an 
individual  instrument  to  be  less  than  a  few  percent,  equation  (17) 
may  be  simplified  as 


2 


(IS) 


It  should  be  explained  that  when  the  polar  difference  is  out¬ 
side  the  critical  region  for  n  =  2,  if  data  closer  to  the  theoretics 
standard  is  selected,  then  the  system  will  fail  only  when  both  data 
are  ineffective.  This  is  an  exception  to  equation  (13). 


TABLE  3 


a  12  3*50  7  8  9  10 


Q, 


O' 


9 


l.Sg*  1 2.5 0*  5 q* 


If  the  rate  of  failure  of  a  single  instrument  is  1% ,  then  with 

r  it  +  i  i 

n  instruments,  the  rate  of  failure  will  be  lowered  by  2L  2~J— 1 
order  of  magnitude. 

3.  The  accuracy  in  averaging  the  normal  data  after  polar  difference 
checking . 

Element  failure  may  be  divided  into  two  types.  Cne  is  sudden 
failure,  meaning  that  a  fast  rotation  is  caused  by  the  failure  of 
one  part.  The  other  is  property  failure,  meaning  that  error  may 
become  larger  than  the  critical  range  W  due  to  the  accumulation  of 
random  interference  with  the  same  sign.  The  former  is  discussed 
under  reliability  in  equation  (17),  while  the  latter  is  discussed 
under  accuracy  analysis,  derived  as  follows. 

Let  n  numbers  be  measured.  Two  cases  will  be  differentiated. 

One  case  is  that  polar  difference  checking  is  normal.  The  probabi¬ 
lity  of  occurrence  of  this  case  is  1-a.  After  arithmetic  averaging, 

the  mean  square  deviation  is  o,  :  the  other  case  is  that  Dolar 

1-a 

difference  checking  has  not  passed.  The  probability  of  occurrence 

for  this  case  is  a.  If  the  disqualified  data  are  not  discarded, 

the  mean  square  deviation  after  averaging  is  a  . 

3. 

The  above  two  cases  combine  to  be  ordinary  arithmetic  averaging, 
namely , 

—a)  +  ola  (19) 

Apparently 

°m>  vT><7—  (2G) 

If  polar  difference  checking  has  failed,  after  discarding  the 
disqualified  data,  the  mean  square  deviation  after  averaging  is 
Let  us  discuss  the  worst  condition  for  a  .  The  worst  situation  is 
to  select  the  median  of  the  sequential  data  (i.e.,  arranging  the 
data  sequentially  from  small  to  large,  take  the  middle  value  for  odd 
number  of  values,  and  the  average  of  the  middle  two  values  for  even 
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number  of  data).  In  reference  [4],  it  is  proved  that  n  <  10  the 

worst  situation  .a,  =  l.2-?~  .  From  [6]  cage  421  example  6-2,  it  is 

v  n 

also  proved  that  the  mathematical  expectation  and  square  deviation 
of  a  simple  median  approach  is  respectively  m  and  .  From 

this  one  obtains 


1.1 
V  2 


=  1.25 -C 

N  FI 


(21) 


Hence,  the  worst  combined  mean  square  deviation  a ^  is 
a,  =  Va‘—(1-a)  +l-25!-^-a  ^  ^-[1-a- 1 .25^a] 

If  we  let  a  =  0.1,  then 

0/<l  .025  -%* 

V'  FI 


(22) 
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From  this  it  may  be  concluded  that  polar  difference  checking 
has  greatly  enhanced  reliability,  with  accuracy  approaching  that  of 
arithmetic  averaging,  i.e.,  by  -L,  times  with  the  worst  case  not 
to  exceed  1  • 025 

IV.  Exanoles 


Example  1.  Method  with  three  sets  of  duplications 

Method  1.  The  middle  value  is  selected  directly  without  polar 
difrerence  checking.  For  instance,  integrate  the  vehicle  accelera¬ 
tion  with  three  acceleration  integrators.  When  the  present  velocity 
is  reached,  an  engine  shut-off  command  is  issued.  As  there  exist 
random  errors  in  all  three  integrators,  the  issuing  of  the  engine 
shut-off  commands  will  not  be  at  the  same  time.  Only  when  the  second 
command  has  been  received  will  the  engine  be  actually  shut  down.  Here 
the  computational  setup  may  be  saved,  since  polar  difference  checking 
is  not  necessary. 

Method  2.  Arrange  the  data  measured  In  real  time  from  small  to 


large  as  x1,x2,x^. 

1 

take  x~; 

and  when  < 

c 

j 

take  x - 
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Method  3.  Select  JC=  j(*1-r *,  +  *,)  without  checking, 

methods  are  compared  as  follows: 


The  three 


Y  « 

1 

* 

[T  \ 

i  *  *  m 

y 

r 

I 

1 _ 

0.577 O 

o.tnta 

0.577® 

Q> 

1 

1  4 

_ 1 _ 

*4 

l.Sg‘ 

1.54* 

1 — single  meter;  2 — ordinary  arithmetic  averaging;  3 — median  value; 
4 — polar  difference  checking;  5 — mean  square  difference  a  ; 

6 — failure  rate  Q  ^ 

P 


Example  2.  With  four  duplication  sets,  the  accuracy  with  averaging 

after  polar  difference  checking  is  enhanced  one  fold  comparing  to 

2 

that  of  a  single  meter.  The  failure  rate  is  3q  .  It  should  be 
pointed  out  that  when  there  exist  unknown  regular  errors  for  the 
single  meter,  there  are  two  factors  m  and  o  in  measuring  the  accuracy 
of  the  single  meter.  Now  the  probability  for  satisfying  the  condi¬ 
tion  that  the  maximum  error  is  to  be  within  the  range  Zmax  is  99. 3"- 
If  the  regular  error  m  =  0,  then  the  single  meter  accuracy  is 

required  to  be  a=Z  /2.7.  When  m  /  0,  what  should  be  the  single 

m3  x 

meter  accuracy  a 2  so  that  the  requirement  that  the  probability 

range  be  within  the  allowed  limits  may  be  satisfied?  Let  m  =  rax,  41 

ax  =  Ba,  the  following  equation  should  be  satisfied 


'  1  -<*-y0<T  )'n<0oy 

J 2n  Pa e 

^0 


dx  =  99.3% 


which  may  be  solved  to  get 


V 

0 

os 

0.1 

0 

1 

0.011 

0.124 

With  four  identical  elements,  half  of  it  installed  normally  and 
the  other  half  reversed  and  backward,  the  regular  error  is  basically 
eliminated.  With  polar  difference  checking  average,  the  random  error 
is  halved,  so  that  the  accuracy  is  enhanced  to 
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0 


t.s 


accuracy 

enhancement 


o  t 


0  456 


0.412 


In  space  vehicle  design,  the  guidance  error  caused  by  the  error 
in  the  longitudinal  axis  acceleration  meter  is  one  order  of  magni¬ 
tude  larger  than  in  other  directions.  If  we  install  four  duplicate 
meters  in  this  direction,  the  accuracy  of  the  whole  system  will  be 
equivalently  enhanced  one  fold,  and  the  accuracy  is  also  greatly 
enhanced . 


Example  3-  Realistic  duplication  method  of  strap-down  guidance  sys¬ 
tem.  Duplicate  signals  for  pitch,  yaw  and  roll  angles  may  be  mea¬ 
sured  with  three  gyros  with  two  degrees  of  freedom.  Duplicate  signals 
for  longitudinal  acceleration  may  be  measured  with  two  longitudinal 
acceleration  meters.  A  transverse  acceleration  meter  is  secondarily 
important  and  will  not  be  duplicated.  When  the  duplicate  signal 
polar  difference  exceeds  the  critical  value,  the  data  closer  to  the 
theoretical  standard  will  be  selected.  When  the  polar  difference  is 
less  than  or  equal  to  the  critical  value,  the  average  of  ;he  two  values 
is  selected.  In  this  way,  although  one  acceleration  meter  and  one 
gyro  with  two  degrees  of  freedom  are  added,  the  guidance  error  is 
reduced  by  times,  and  the  failure  rate  by  two  orders  of  magni¬ 

tude.  When  this  kind  of  setup  is  installed  in  vehicles  with  a  fast 
engine  propeller,  there  will  be  no  need  to  change  elements  if  it 
should  be  discovered  before  launching  that  some  inertial  element  may 
be  faulty  so  that  time  will  not  be  lest. 
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THE  RESEARCH  OF  THE  PROJECT  FOR  ADJUSTING  THE 
THRUST  OF  SOLID  ROCKET  PROPULSOR 

Xu  NVencan  ■!*,?  •- 

...  .  (  Nj  C J  O  ' , 

Abstract  '  0 

A  conception  of  adjustable  solid  propellant  rocket  motor  was  proposed.  It  is  going  to  be 
achieved  through  the  use  of  propellants  having  negative  pressure  exponent  n.  In  this  paper, 
the  feasibility  of  the  technique  and  characters  of  adjustment  were  discussed  Other  advantages 
of  the  propellent  with  a  negative  pressure  exponent  were  also  presented. 


SYMBOL  TABLE 

A.  combustion  surface  area 

b 

A*.  jet  throat  area 

a  combustion  rate  coeffi¬ 

cient 

C*  characteristic  speed 

C_  thrust  coefficient 

C.(i=l,2)  constants 

F  thrust 

?  maximum  thrust 

max 

?  given  thrust 

g 

k  specific  heat  ratio  of 

fuel  vapor 

X, (1*1,2)  constant  coefficients 
1  (including  those  with 

superscript 

n  combustion  pressure 

index 


P  combustion  chamber  pressure 

c 

r  combustion  rate 

t  time 

u(t)  unit  step  function 

V  combustion  chamber  volume 

c 

a  (=X/pc )  propellant  thermal 

diffusion  rate 

T  function  of  k 

e(»)  static  difference 

X±(i=l,2)  constant  coefficients 

p  fuel  vapor  density 

c 

Pp  propellant  density 

o  superharmonics 

To  o c  transition  time 
0  .  Up 


^3 


Introduction 


Much  work  has  been  done  both  internationally  and  domestically  on 
the  thrust  control  of  solid  propellant  rocket  engines,  resulting  in 
the  development  of  such  direction  control  systems  and  thrust  termin¬ 
ation  mechanisms  based  on  various  action  principles  as  "inclined 
flow  ring",  "rotating  Jet",  "vibrating  jet",  "soft  jet",  "liquid 
float  jet"  and  "secondary  injection"  of  liquid  of  fuel  vapor,  etc. 
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However,  there  is  as  yet  no  good  and  practical  method  for  regulat¬ 
ing  the  strength  of  thrust.  Although  there  have  been  some  proposals 
such  as  "layered  engine"  to  inject  certain  liquids  into  the  combus¬ 
tion  chamber  and  some  preliminary  attempts  in  applying  methods  of 
installing  a  regulating  cone  in  the  jet  or  making  use  of  propellant 
with  n  greater  than  0  and  close  to  1  with  a  mechanical  rotating 
valve  based  on  propellants  with  positive  pressure  coefficients  [1] 

[2],  yet  little  prospect  can  be  seen  due  to  the  complexity  in  the 
structure  and  the  bad  weight  and  service  properties. 

In  this  paper,  we  shall  attempt  to  investigate  the  possibility 
of  solving  this  problem  by  using  new  propellants  with  pressure  coeffi¬ 
cient  less  than  0,  but  absolute  value  greater  than  1  (for  convenience 
of  discussion,  we  shall  call  them  by  the  name  "high  negative  pressure" 
propellant ) . 


The  feasibility  and  advantage  of  implementing  thrust  control 
"high  negative  Dressure  propellant" 


Why  do  we  want  to  study  this  kind  of  propellant?  Will  this 
kind  of  propellant  provide  the  feasibility  of  implementing  thrust 
regulation?  What  are  the  advantages?  Apparently,  these  questions 
should  be  answered  first.  As  is  well  known,  when  the  combustion  rate 
obeys  the  equation  r  =  ap„,  and  the  sound  speed  flow  exists  in  the 
jet  throat,  the  balanced  pressure  and  the  corresponding  thrust  of 
the  solid  fuel  rocket  engine  combustion  chamber  may  be  expressed 
with  the  two  equations  below: 

P.  =  (  (1) 

F  =  C,P^i,  (2) 

From  the  equations  above,  it  can  be  seen  that  once  the  propell¬ 
ant  has  been  selected,  the  only  effective  way  to  regulate  the  thrust 
is  to  vary  the  parameters  Afe  and  At .  Of  the  two,  the  parameter  that 
allows  one  to  randomly  vary  it  by  simple  means  (i.e.,  the  kind  of 
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parameter  that  is  varied  not  according  to  a  preset  program,  but 
according  to  current  need)  is  basically  the  jet  nozzle  throat  area 
A t .  Hence,  we  shall  vary  equations  (1)  and  (2)  simultaneously  to 
obtain  the  equations  below  as  the  basic  equations  in  this  paper: 


PC  =  K,A,”-' 

(  3  ) 

(4) 

P,  =  K,F''' 

(  5  ) 

where 


K i  —  (P  f  •  C*’Ab-  a)  i-n 

K,  =  C,-Kt 

.  K,~  • 


When  the  propellant  is  already  determined  and  the  combustion 
area  is  a  constant,  it  is  apparently  permissible  to  use  K,  as  a 
constant  in  the  analysis.  Cp  is  also  obviously  a  function  of  A^/A^ 
and  ?  ,  but  in  the  confines  of  our  discussion  here,  it  may  be  taken 
approximately  as  a  constant.  Thus,  we  shall  treat  K2 ,  approxi¬ 
mately  as  constants  below.  In  this  way,  by  differentiating  logar¬ 
ithmically  both  sides  of  the  equations  above,  we  obtain 


dp,  ^ 1 d.A , 

p.  *ii-1  A, 

dA ,  1  \dF 

A,  \  n  )  F 

dp,  ^  1  dF 

P,  *  »  F 


(  7  ) 

it) 
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From  equations  (5)  or  (S),  we  know  that  in  order  to  regulate  the 
thrust,  it  will  be  necessary  to  alter  the  pressure  in  the  combustion 
chamber.  But  when  the  variation  in  pressure  becomes  unacceptable  in 
practice,  the  regulation  of  the  thrust  will  also  become  an  empty 
dream.  Hence,  we  must  try  to  minimize  the  pressure  variation  of 
the  combustion  chamber  during  thrust  regulation  (thus  taking  full 
advantage  of  the  structural  strength  of  the  engine),  and  maximize 
the  absolute  value  of  n.  For  example,  when  the  pressure  index  n  = 

0.75,  the  pressure  must  be  increased  16  times  or  6.3^  times, 
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respectively,  in  order  to  increase  the  thrust  by  8  or  4  times. 
Changes  of  this  magnitude  will  naturally  lower  the  weight  property 
of  the  engine  when  we  want  to  make  sure  that  it  will  operate  under 
the  highest  pressure.  This  is  also  the  point  of  vulnerability  in 
some  of  the  experiments  and  design  reported  in  some  countries. 

But  when  n  =  -2,  to  increase  the  thrust  by  the  same  factor,  the 
pressure  needs  to  be  reduced  only  to  1/3. 83  or  1/2  of  the  original 
pressure.  Obviously,  the  situation  is  greatly  improved  and  the 
combustion  chamber  pressure  variation  significantly  reduced. 


From  the  above  consideration,  it  is  apparent  that  the  larger 
| n |  is,  the  better.  However,  in  order  to  keep  the  engine  operating 
stably,  a  n  value  with  absolute  value  greater  than  1  cannot  be 
positive  because  when  the  flow  rate  of  the  engine  is  equal  to  the 
rate  of  generation  of  the  fuel  vapor,  it  operates  under  equilibrium 
pressure  (o,^)  and  for  this  equilibrium,  like  all  other  equilibrium 
problems,  there  exist  stable,  unstable  and  neutral  states.  If  due 
to  the  effect  of  some  random  factor  the  combustion  chamber  pressure 
deviates  from  the  equilibrium  pressure  Pcb,  then  will  there  be  a 
tendency  to  diminish  this  deviation,  or  a  tendency  to  increase  this 
deviation?  Or  will  equilibrium  be  established  under  the  new  condi¬ 
tions?  These  are  then  the  difference  between  the  three  types  of 
equilibria.  For  the  engine  to  operate  stably,  we  must  have  stable 
equilibrium.  The  mathematical  expression  for  this  analysis  may  be: 
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Hence  in  order  that 


JL-(  -d£-)<o 

dp,  \  it  > 


We  must  let  n  -  1  <  0  46 

i .e . ,  n  <  1 

This  may  be  seen  intuitively  from  Figure  1. 


Jf 


Figure  1 

anslator’s  note:  illegible 


Figure  2 
in  foreign  text. 


Hence,  we  should  use  "high  negative  pressure  index"  propellant 
so  as  to  minimise  the  variation  of  the  combustion  chamber  pressure 
during  thrust  regulation. 


a) 


;  ) 


is  just  the  common  mathematical  expression 

->»•  since  the  former  may  become 


d>r 


9 


At  the  came 
logarithmically , 


time,  differentiating  both  aides  of 
we  find  from  equations  (1),  (2)  that 

=  _J _ r  dPP  dC*  ,  dAt  dA,  d.  I 

1-n  L  P,  C*  '  .4k  .4,  +T-T-J 


equation  (1) 

(f) 


It  can  be  seen  that  fluctuations  In  the  parameters  p,X*.-4*.  A,,  a 
induced  by  some  factor  will  also  cause  the  combustion  chamber  press¬ 
ure  p  and  the  thrust  F  to  fluctuate  and  when  the  relative  fluctua- 
c 

tion  of  the  former  is  fixed,  that  of  the  latter  will  decrease  with 
decreasing  n  value.  In  practice,  among  the  above  parameters  (i.e., 
P,X*.Ak%A1%a  the  combustion  rate  coefficient  a  fluctuates  the  most. 
Based  on  the  demand  of  practical  fuel  rod  production,  it  is  often 
permitted  to  vary  over  a  wide  range.  For  example,  in  the  two  typical 
ways  of  preparing  mixture  propellant,  the  values  of  Aa/a  are  speci¬ 
fied  respectively  to  be  about  0.07  and  0.02  which  is  already  fairly 
difficult  in  the  industry.  But  if  n  =  0.75,  then  together  with  the 
fluctuations  of  ether  parameters,  Apc/pc  may  reach  as  high  as  0.28 
and  over  0.08,  respectively.  However,  if  n  =  -2,  then  Apc/pc  will 
be  respectively  reduced  to  0.023  and  0.007.  The  advantage  is  obvious 
It  can  significantly  decrease  the  fluctuation  in  the  combustion 
chamber  pressure  and  the  spread  of  engine  thrust. 


Besides  the  effect  of  kA^/A^  on  Apc/pc  can  also  be  reduced 
significantly,  thus  permitting  the  design  of  fuel  rods  with  large 
area  ratio  to  increase  the  fuel  Insertion  density  of  the  engine. 


Summarizing  the  above  discussions,  as  shown  in  Figure  2,  by 
using  a  new  "high  negative  pressure  index"  propellant,  it  is  possi¬ 
ble  to  regulate  thrust  within  a  wide  range  without  altering  the 
current  basic  structure  of  the  solid  propellant  rocket  engine.  It 
also  possesses  the  advantage  of  reducing  the  fluctuation  in  combus¬ 
tion  chamber  pressure  and  the  spread  in  engine  thrust  and  increasing 
the  density  of  fuel  rod  insertion. 


Of  course,  in  order  to  implement  the  regulation  of  thrust,  it 
is  also  necessary  to  be  able  to  regulate  over  a  wide  range  the  area 
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of  the  nozzle  throat.  It  can  be  observed  from  equations  (4)  and 
(7)  that  when  n  <  0  and  |n|  >  1,  one  needs  to  vary  the  nozzle  throat 
area  A  over  a  large  range  to  obtain  the  required  thrust  regulation 
range.  In  the  two  examples  cited,  the  thrust  is  increased  eight¬ 
fold  in  both  cases.  For  n  =  0.75,  A  needs  only  be  reduced  by  one- 
half,  while  for  n  =  -2,  A  needs  to  be  increased  by  more  than  22.6 
times  its  original  value.  This  is  Inconceivable  in  the  ordinary 
design  of  the  mechanical  regulation  of  nozzle  throat  area.  But  this 
is  not  a  problem  when  "vortex  flow  valve"  of  the  flow  control  type 
is  used.  As  the  energy  of  the  propellant  is  improved  constantly, 
it  is  only  through  the  flow  control  technology  that  the  nozzle 
throat  area  (hence  the  flow  volume)  may  be  regulated  under  the  bad 
'working  condition  of  the  nozzle  throat.  Hence,  theoretically  this 
need  is  not  an  impediment  to  the  implenetation  of  thrust  regulation 
in  our  design.  Cn  the  contrary,  it  actually  provides  the  feasibility 
that,  as  one  part  of  the  automatic  regulation  system,  it  will  achieve 
the  accuracy  and  static  characteristics  desired. 

Ill .  Can  "high  negative  pressure  index"  propellant  be  obtained? 

Limited  by  the  study,  we  are  not  yet  able  to  give  a  quantitative 
analysis  based  on  combustion  mechanism.  We  shall  only  mention  the 
following  points  as  a  basis  in  our  discussion: 

1.  Propellants  with  n  <  0  already  exist  not  only  theoretically,  but 
also  in  practice.  In  one  preparation,  a  pressure  range  of  50-70 
kg/cm  and  pressure  index  n  =  -0.825  has  been  achieved. 

In  the  reports  on  double-base  platform  propellant  and  mixed 
material  propellant  overseas,  the  so-called  "mesa  effect  [3,5]  with 
n  <  0  has  appeared.  Recently,  there  have  been  reports  on  mixed 
material  propellant  with  n  <  0  and  even  n  =  -4.0  [4].  This  demon¬ 
strates  that  it  is  feasible  and  practical  to  obtain  "high  negative 
pressure  index"  propellants. 
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2.  In  many  foreign  references,  [1,4,5,61,  the  mechanism  of  the 
"mesa  effect"  in  the  combustion  of  mixed  material  solid  propellant 
has  been  discussed. 

(1)  From  the  microfilm  of  the  combustion  surface  after  exper¬ 
imental  ignition  shut-off,  it  can  be  observed  that  at  low  pressure 
(less  than  30  kg/cm),  the  surface  of  ammonium  chlorate  (Ap)  crystals 
is  convex,  since  the  thermal  dissociation  of  the  adhesive  under 
this  condition  is  faster  than  that  of  Ap,  while  at  high  pressures 
(greater  than  30  kg/cm),  the  surface  has  concave  pits,  indicating 
that  the  Ap  crystals  are  consumed  faster  [1]. 

(2)  p-amino  ester  (pU)  is  different  from  the  adhesive  carboxyl 
p-butadiene  (CTPB)  in  that  it  already  becomes  very  fluidized  before 
the  temperature  for  its  rapid  dissociation  is  reached,  while  the 
latter  is  a  highly  adhesive  foamy  fluid.  Under  high  pressure,  local 
quenching  [1,4,5,61  will  occur  when  an  easily  fluidized  adhesive 
flows  over  and  covers  the  surface  of  Ap  crystals  that  dips  below 
the  combustion  surface. 

(3)  It  has  been  observed  with  the  movie  camera  method  that 
irregular  pulsed  combustion  with  partial  quenching  has  appeared  in 
platform  propellants.  The  area  of  this  local  quenching  as  a  frac¬ 
tion  of  the  total  area  determines  the  degree  of  decrease  of  the  com¬ 
bustion  rate  averaged  over  the  whole  combustion  surface.  The  depen¬ 
dence  of  this  fraction  on  the  pressure  is  what  causes  the  decrease 
of  the  pressure  index  n  until  the  appearance  of  the  "mesa"  effect 
with  n  <  0  [1,4,51. 

(4)  When  salt  that  melts  easily  such  as  (NH^KSO^  is  added 
to  the  propellant  as  combustion  speed  regulator,  the  local  quench¬ 
ing  phenomenon  can  also  be  induced.  Or  the  dissociation  of  NK^CIO^ 
may  be  suppressed  [41  due  to  the  NH*  ions  and  a  deep  thermal  poten¬ 
tial  hill. 
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3.  In  [4],  a  physical  picture  verified  by  some  expeirmental  obser¬ 
vations  is  proposed  for  the  combustion  process  of  this  kind  of  pro¬ 
pellant  . 

An  engine  experiment  using  80  lb.  of  propellant  with  n  =  -2.5 
is  also  reported. 


IV .  Preliminary  analysis  of  some  regulatory  properties  4 8 

For  the  overall  cycle  of  thrust  auto-regulating  systems,  the 
regulation  of  the  engine  is  only  one  link.  However,  its  regulatory 
properties  obviously  should  be  analyzed. 


For  convenience,  let  us  assume  that  the  combustion  rate  still 
obeys  the  law  r  =  ap((  during  the  regulatory  process,  that  at  all 
times  sound  velocity  flow  exists  at  the  nozzle  throat,  and  that  the 
pressure  in  the  combustion  chamber  is  the  sane  everywhere  and  is  also 
equal  to  the  pressure  at  the  nozzle  inlet.  Thus,  we  get  from  the 
continuity  equation  of  the  fuel  vapor: 


dP, 

dt 


(10) 


Again  by 
constant 


and 


differentiating  equation  (2)  logarithmically  with  CP  a 
we  get 

i IF  _  a  p,  dA, 

F  p,  ~A,~ 


1  dF  _  1  cl  P  1  J.l, 

F  dt  p.  dt  A,  at 


<  ID 


5y  solving  simultaneously  equations  (2),  (10)  and  (11),  we  can 
obtain  the  dynamic  equation  for  this  part  as 


_1  dF 
F~  dt 


—  F,A, 


l_  _dA, 
A,  dt 


(12) 


where 
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r*c ** 

A.=  (P,-Pe).ifaCy 

k-  A’C* 

A  ‘“T. 

When  the  input  is  a  unit  step  function,  i.e.,  A^  =  A^  +u(t) 

(Afa  is  the  nozzle  throat  area  before  the  unit  step  functiBn  action) 
apparently  dA^/dt  =  0  and  u(t)  =  1  for  t  >  0.  By  differentiating 
equation  (12),  we  get 

^-  +  K;F.K\F- 

where 

y.  K. 


The  differential  equation  obtained  for  this  unit  transition 
process,  equation  (13),  is  a  typical  Bernoulli  equation.  Hence, 
its  solution  may  be  obtained  through  the  method  of  separation  of 
variables.  When  n  /  1.0, 

~TT-^nTK r  1 n ( K '  ~  K  ,F  1 " " )  M  c ■ 

From  the  initial  conditions:  t  =  0+,  F  =  F(0),  we  get 

C,=  (l-n)Ai  'o  K'~K^F(0)j  •"} 

Therefore , 


where 


t  = 


_ 1 _  ]n  A  i  —  K 

(l-n)h’j  KV-K'tFr-' 


i 

7r(0  =  A,{l  —  A,expC(n  —  *“• 


a:  =  i-(a;/a;)[F(o)]'" 


Also 

and 


since  (n-1)  <  0,  then  from  (19)  we  get,  when  t  - >°°,  F(ot)  =  X 

A  - i-f ^r(c°)T'1  .  Hence  eauation  (14)  may  also  be  written  as 

1  F(0)  J 


F(O**F(0)j 


(F( 0)  v" 

\F(  oT)/ 


i  expC(fi-l)  /f^]}itrr 


C/s') 
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1.  Transition  period:  by  definition 
we  get 


lF(r».t>)-F(ce)i 

F(o=) 


=  0.05 


.  _  1  ,  1-(1±0.05)'-" 

0,01  (w— l)Ki  .  .  F(0)  r* 

- 


C/O 


where  the  +  and  -  signs  correspond  respectively  to  the  situations 
of  0  <  n  <  1  and  n  <  0. 

When  the  relative  regulation  is  assumed  to  be 

=  =  10<$ 

F(oo) 

for  the  cases  n  =  0.75  and  n  =  -2,  tq  is  respectively  about 
1.43/K£  and  0.213/K1,  which  means  that  the  transition  period  of  the 
latter  is  only  1/7  of  that  of  the  former. 


2.  Super-regulation:  From  equation  (14),  we  get  Fmax  =  X-,  •  There¬ 
fore,  by  definition  the  super-regulation  is 

o  _  .  i oo *  - o . 

F(o°) 


Static  difference 


-4k 


and 

we  car.  cbtain^»  =  j  i-»’*F(oc) 


he  given  value  =  ^ r/fr[^r,TU(i)] 
i]t=T. 

Therefore,  by  definition,  the 


static  difference  .ioo»»o  .  (In  the  above  F  is 

/  F,  eg 

the  static  equilibrium  pressure  corresponding  to  a  given  thrust). 


When  the  pressure  changes  rapidly,  the  combustion  rate  is  no 

longer  only  a  function  of  the  pressure,  but  is  also  a  function  of 

the  rate  of  change  cf  the  pressure,  e.g.,  r*«ap;  (\+ — ®2__  dpt  \ 

'  a'p/*  dt  r 

(If  we  assume  that  the  empirical  equation  can  still  be  used  for  n  < 

0,  we  can  obtain  by  the  same  token  the  corresponding  dynamic  equation 
for  unit  step  process: 

dF  _  K7F*"-  -KtF"' 
dt  K"F'"  —  ic%m 


When  n  /  1.0,  its  solution  is 


<  =  ■ 


(1  -n)K\  ln(A'  K'F"">  /*1"!  ) 


Experiments  indicate  [7]  that  the  combustion  rate  pressure  index 
n_  for  instantaneous  pressure  change  will  still  be  greater  than  the 

Cl 

pressure  index  n  during  slow  change.  Hence,  the  assumption  above 
that  the  combustion  rate  during  the  regulatory  process  agress  with 
t  =  apj  is  obviously  going  to  be  quite  different  from  reality.  This 
difference  may  be  seen  from  the  empirical  law  shown  in  Figure  2  of 
reference  [4].  Therefore,  the  above  analysis  is  only  preliminary. 
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